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ABSTRACT 

We present the Planck Sky Model (PSM), a parametric model for the generation of all-sky, few arcminute resolution maps of sky emission at 
submillimetre to centimetre wavelengths, in both intensity and polarisation. Several options are implemented to model the cosmic microwave 
background. Galactic diffuse emission (synchrotron, free-free, thermal and spinning dust, CO lines). Galactic Hii regions, extragalactic radio 
sources, dusty galaxies, and thermal and kinetic Sunyaev-Zeldovich signals from clusters of galaxies. Each component is simulated by means 
of educated interpolations/extrapolations of data sets available at the time of the launch of the Planck mission, complemented by state-of-the-art 
models of the emission. Distinctive features of the simulations are: spatially varying spectral properties of synchrotron and dust; different spectral 
parameters for each point source; modeling of the clustering properties of extragalactic sources and of the power spectrum of fluctuations in the 
cosmic infrared background. The PSM enables the production of random realizations of the sky emission, constrained to match observational data 
within their uncertainties, and is implemented in a software package that is regularly updated with incoming information from observations. The 
model is expected to serve as a useful tool for optimizing planned microwave and sub-millimetre surveys and to test data processing and analysis 
pipelines. It is, in particular, used for the development and validation of data analysis pipelines within the Planck collaboration. A version of the 
software that can be used for simulating the observations for a variety of experiments is made available on a dedicated website. 

Key words. Cosmology: cosmic background radiation - Interstellar medium (ISM), nebulae - Galaxies: clusters: general - Galaxies: general - 
Infrared: diffuse background 



1. Introduction 

The cosmic microwave background (CMB), relic radiation from 
the hot Big Bang, carries an image of the Universe at an age 
of about 380,000 years. CMB anisotropies reflect the inhomo- 
geneities in the early Universe, and are thus an observable of 
great importance for constraining the parameters of the Big Bang 
model. 

For this reason, a large number of experiments dedicated to 
CMB anisotropy detection and characterisation have observed 
the sky at wavelengths ranging from the sub-millimetre to a few 
centimetres (frequencies from a few to a few hundred GHz). 
Stimulated by CMB science, astronomical observations in this 
wavelength range have enriched many fields of scientific inves- 
tigation, ranging from the understanding of emission from the 
Galactic interstellar medium (ISM), to the study of a large popu- 
lation of extragalactic objects, including radio galaxies, infrared 
galaxies, and clusters of galaxies. 

The WMAP sateflite, launched by NASA in June 2001, has 
surveyed the sky in five frequency bands, ranging from 22 to 
94 GHz ( [Bennett et al.|[2003 ). The data from this mission have 
provided an all-sky map of CMB temperature, and an unprece- 
dented data set for the study of sky emission at millimetre wave- 
lengths -both in total intensity and in polarisation ( Jarosik et al. 



2011| ). The Planck missioiQ launched by ESA in May 2009, has 
started to provide even better observations of the sky in nine fre- 
quency bands, ranging from 30 to 850 GHz ( Tauber et al. |2010t 
[ Planck Collaboration I[ |201 These observations, which will 
be released to the scientific community early 2013, are expected 
to become the new reference in this frequency range. 

It has long been recognised that planning future observations 
and developing methods to analyse CMB data sets both require 
realistic models and simulations of the sky emission and of its 
observations (Gawiser et al. 1998; Bouchet & Gispert 1999 



Giardino et al.„2002[]de Oliveira-Costa et al.,,2008„Sehgal et al. 



|2007 2010| l.For the preparation of the analysis of CMB data sets 
such as those of the Planck mission, it has proven useful to put 
together a model of sky emission, the Planck Sky Model (PSM), 
which can be used to predict or simulate sky emission for various 
assumptions and various parameter sets. The pre-launch version 
of this model, based on publicly available data sets existing be- 



fore Planck observations (see Section 2.3 1, is used in particular 



as a simulation tool in the context of Planck. An update of the 
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Planck (http://www.esa.int/Planck) is a project of the European 
Space Agency - ESA - with instruments provided by two scientific 
Consortia funded by ESA member states (in particular the lead coun- 
tries: France and Italy) with contributions from NASA (USA), and tele- 
scope reflectors provided in a collaboration between ESA and a scien- 
tific Consortium led and funded by Denmark. 
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model is planned on the basis of observations obtained by the 
Planck mission. 

This paper summarises the outcome of this preparatory mod- 
elling and simulation effort. Software and simulations are made 
available to the scientific community on a dedicated websit^ 
The current version of the Planck Sky Model (PSM), described 
in the present paper, is vl.8. We plan regularly to update mod- 
els and software tools on the basis of upcoming observations, 
and when additional sophistication becomes useful for upcom- 
ing experiments and scientific investigations. The PSM devel- 
opers welcome suggestions for improving the models or imple- 
menting alternate options. 

Although nothing prevents running the PSM to simulate sky 
emission at any frequency, our sky model is actually valid, at 
present, only for frequencies ranging from about 3 GHz to about 
3THz, i.e., wavelengths ranging between 100 //m and 10 cm. 
At lower frequencies, Faraday rotation is important and our po- 
larised maps are not representative of the sky emission, as this 
effect is presently neglected. Intensity maps, however, should 
be reasonably accurate down to about 500 MHz. At frequencies 
higher than 3 THz, the complexity of the dust emission is not 
properly taken into account by the present model, for which the 
representation of dust emission is based on spectral fits valid be- 
low 3 THz. 

The paper is organised as follows: in Section 2 we review 
the basic features of the model as well as the architecture of the 
package; in Sections 3, 4 and 5 we describe the simulation and 
implementation of CMB anisotropics, diffuse Galactic and com- 
pact object emissions; and finally, in Section 6, we summarise 
the status of our sky model at the time of this publication and 
draw some conclusions. 



2. The Model 

2.1. Rationale and examples of application 

The development of a sky model and sky-emission simulations 
has several complementary objectives. First, a prediction of sig- 
nals is useful for the optimisation of planned instruments (in 
terms of resolution, number and characteristics of frequency 
channels, acceptable noise level, and in general for making the 
necessary compromises between instrumental performance and 
costs or feasibility). Secondly, an emission model is necessary 
to simulate representative 'observed' data sets, i.e., simulate ob- 
servations of the modelled sky, which can then be used as test 
cases for developing dedicated data analysis software and for 
testing data analysis pipelines. Finally, such a model is a tool for 
easy comparison of any observed data set (at the time it becomes 
available) with what is expected based on previous empirical 
or theoretical understanding. Because of these different objec- 
tives, two different philosophies co-exist in the development of 
the model. 

For some applications, we address the problem of making an 
accurate prediction of sky emission, considering existing obser- 
vations and present knowledge about the emission mechanisms. 
The objective is then to give the best possible representation of 
our sky as it has been observed. Such a prediction can be used as 
a basis for calibration of future data against a common model, for 
selecting sky areas to be observed or for cleaning observations 
from expected contamination by a particular component. A pre- 
diction of sky emission is strongly data-driven; when only upper 
limits are available on the basis of existing observations (e.g.. 



CMB B-modes, cluster peculiar velocities, . . . ), we assume our 
knowledge to be non-existent, and a prediction is not possible. 

For other applications, such as the development of data pro- 
cessing and analysis pipelines, or comparison of methods such 
as in Leach et al. ( 2008| , we need a realistic, statistically rep- 



http://www.apc.univ-paris7.fr/~delabrou/PSM/psm.html 



resentative, sky simulation, with the appropriate complexity and 
with some level of randomness, compatible with current uncer- 
tainties in the observational data and their interpretation in terms 
of a model. 

Our model is designed to implement both sky emission pre- 
dictions and simulations. For practical uses, the two are com- 
bined to achieve a realisation of the sky, at the desired com- 
plexity and resolution, in the context of a given experiment or 
general study. For each of the components, several options exist 
for the modelling, this feature allowing the user to investigate 
the impact of theoretical uncertainties in our interpretation and 
parametrisation of sky emission. The PSM can also be used to 
generate model data sets similar to those which could be pro- 
duced by an instrument. Basic 'sky observation' tools allow for 
the integration of the sky emission in frequency bands, smooth- 
ing with instrumental beams, and observing along scans. A com- 
parison of 'observed' sky predictions with actual data being col- 
lected by upcoming instruments allows for real-time assessment 
of the proper operation of experiments. 

The sky simulation tool implemented with the PSM was 
originally developed for investigating component separation 
methods, building on early work by |Bouchet & Gis pert (1999]l 
for the Phase A study of the Planck space mission. Sky obser- 
vations targeting the measurement of CMB intensity and polar- 
isation anisotropics in fact contain a superposition of emission 
from several other astrophysical sources: the diffuse interstellar 
medium of our Galaxy; radiogalaxies and infrared galaxies; the 
Sunyaev Zeldovich (SZ) effect from the hot intra-cluster electron 
gas in massive clusters of galaxies. Separating these components 
into contributions from distinct processes is important for the 
interpretation of the observations. In the past 15 years, compo- 
nent separation methods developed by several authors have been 
making use of simulated sky maps to assess the performance of 
the proposed approaches, or to validate them before application 
to real data sets. The large number of corresponding publica- 
tions testifies of the importance of this field of investigation (see, 
e.g.|Baccigalupi et aL] [2000; Sno ussi et al.||2002^|Delabrouille, 
Cardoso, & Patanchon, 2003; Martinez-Gonzalez et al.' '2003' 
Bedini et al. 2005 ; Stolyarov et al. , 20 05 ; Patanchon et al. , 2005 
Bonaldi et al. 2006; Eriks en et al.| |2 006; StivoH et al. 2006 
Bob in et al., ,2008 ; Leach et al.[ 12008, ,Stompor et al.) ,2009 
N0rg aard-Nielsen & Hebert | 2009 [Efstathiou, Gratton, & Pa 
2009; Betoule et al. 2009; Dick, Remazeilles, & Delab rouille 
2010; Remazeilles, Delabrouille, & Cardoso, ,201 la. , to mention 
only a fraction). 

It is likely that astrophysical confusion (i.e., contamination 
by other astrophysical emission), rather than by instrumental 
noise, will set the limit on the performance of upcoming CMB 
observations, such as those of the Planck mission. For this rea- 
son, component separation is a very active topic of research, for 
which accurate predictions or realistic simulations of the sky 
emission, such as those provided by the present model, are nec- 
essary. 

Sky emission maps from various astrophysical components 
in the PSM can be generated for different models, all of which 
are regularly updated. For each model, the sky emission depends 
on a set of parameters. The generation of sky maps within the 
PSM allows for comparison of such variants with actual obser- 
vations, and hence for discrimination between models, or their 
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parameters. The CosmoMC tool developed by [Lewis & Bridle] 
( [200 2 ) is an approach for constraining cosmological parameters 
from CMB observations. In the same spirit, the PSM could be 
used in Monte-Carlo simulations to constrain parameters which 
govern the physics of other sky emission mechanisms. 

2.2. Astrophysical components 

The total emission of the sky in the 3 GHz - 3 THz frequency 
range is modelled as the sum of emission from different pro- 
cesses, which we classify for convenience into three large cate- 
gories: 

1. cosmic microwave anisotropics (including the dipole); 

2. Galactic diffuse emission; 

3. emission from compact objects (external galaxies, clusters of 
galaxies, and Galactic compact sources). 

Each identified process with a given emission law is consid- 
ered as an independent astrophysical component. The current 
software implements the CMB dipole, CMB anisotropics, syn- 
chrotron, free-free, thermal dust, spinning dust, CO molecular 
lines, thermal SZ effect, kinetic SZ effect, radio sources, infrared 
sources, far infrared background and ultracompact Hii regions. 
Solar-system emission from the planets, their satellites, from a 
large number of small objects (asteroids), and from dust particles 
and grains in the ecliptic plane (which generate zodiacal light), 
are not modelled in the current version. 

The emission of all components is represented using para- 
metric forms. Some components are modelled with a fixed emis- 
sion law which does not depend on the sky location. In the 
current implementation, this is the case for the CMB, the non- 
relativistic SZ effects (thermal and kinetic), free-free and spin- 
ning dust emission. Synchrotron and thermal dust are modelled 
with emission laws which vary over the sky. Each point source 
is also given its own emission law. 

A summary of the astrophysical components included in the 
model is given in Table [T| 

2.3. Data sets used in the model 



The models of sky emission implemented in the PSM rely on 
maps and object catalogues obtained from observations made 
with different instruments. The main data sets serving as a basis 
for this are full-sky (or near full-sky) observations of the Galactic 
diffuse emission at various frequencies, and catalogues of ob- 
served compact objects. The selection of the data used in the 
model is based on a practical compromise between availability 
and usefulness for the present implementation, and will evolve 
with future releases. 

We currently use: maps from IRAS, DIRBE and WMAP, as 
detailed below; the 408 MHz survey of |Haslam et al. 
maps 



(1982i; 



of CO emission from Dame, Hartmann, & Thaddeus 
( |2001| l; and the Ha maps of Finkbeiner (2003 1. Most of these 
data sets, in HEALPix format (Gorski et al. 2005 1 ), have been 
downloaded from the LAMBDA web sit^ 

The modelling of compact sources makes use of several ra- 
dio surveys (listed in Section 5.2 Table [2]l , of the IRAS cata- 
logue of infrared sources ( jBeichman et al.| 1988 1, and of cat- 
alogues of galaxy clusters observed with ROSAT (compilation 
201 l[l, as w ell as with the Sloan Digital 



from PifFaretti et al. 



of each of the above data sets are postponed to the description of 
individual components below. 

2.4. Sky prediction 

Rather than simply interpolating existing observations, sky- 
emission predictions with the PSM implement recipes for com- 
bining available data to produce a set of modelled astrophysi- 
cal components. Model predictions hence decompose observa- 
tions into components, each of which is described according to 
a parametric model. Such predictions are deterministic, always 
returning the same products once global parameters (frequency, 
pixelisation scheme, coordinate system, model used to interpret 
the data) are set. 

The accuracy of predictions is limited by two factors: the 
quality of sky observations (resolution, noise levels and system- 
atic effects); and the appropriateness of the model assumed to 
make the predictions. 

Predicted maps of emission at a given frequency are de- 
signed to best match the real sky emission as would be observed 
at a target resolution specified by the user The effective resolu- 
tion of the map, however, depends on the available data used to 
make the prediction. Prediction makes no attempt to extrapolate 
observations to smaller scales, whereas simulation does. 

Predicted synchrotron maps are at present determined by 
408 MHz observations and by WMAP observations in the low- 
est frequency channels, and hence have a resolution of about 1°. 
Predicted dust maps, based on Model 7 of Finkbeiner, Davis, 
]& Schlegelj (fT999|, have a resolution of about 9' (except in the 
4 percent of the sky not covered by IRAS, where the resolution 
is that of DIRBE, i.e., about 40')- The predicted CMB is ob- 
tained from WMAP 5-year data, and has a resolution of order 
28', after applying a Wiener filter to minimise the integrated er- 
ror in the reconstructed CMB map (see Section 3.2.2 1. A pre- 



dicted SZ effect is obtained from known ROSAT and SDSS clus- 
ters, for which a model of SZ emission is produced from the 
universal cluster pressure profile, as presented by 



Arnaud et al. 

( 2010| l. The scaling relations for estimating the SZ flux are nor- 



malised on the basis of WMAP and Planck data (pVIehn et al. 



20nj [Piaiick CoUaboration X| [SoTT] [Planck Co llaborat ion XI 
201 II). This prediction contains more than 15,000 clusters (see 
details in Section [5.1.1| l. 

2.5. Simulations 



Simulations generated with our model are random (or partially 
random) realisations. Each time a simulation is generated, a new 
sky is produced. It is possible however to re-generate exactly 
the same sky by rerunning the same version of the PSM, with 
the same sky model options and the same simulation seed (for 
random number generation). 

For given cosmological parameters, the power spectrum of 
the CMB can be predicted using, e.g., the CAMB softwar^ 
(Lewis, Challinor, & Lasenby 2000 1, so that a plausible CMB 



Sky Survey (SDSS, |Koester et al.||2007| l. Details about the use 



sky can be simulated by drawing at random the flf,„s of a sim- 
ulated CMB, according to normal distributions wi th variance 
C(. An interface to the CLASSj^ software (Bias, Lesgourgues, 
|& Tram[|201 l| l is also implemented in the PSM. 

Specific realisations of sky emission can be produced for 
other components as well. A model of cluster counts as a func- 
tion of mass and redshift, dN{M,z)/dMdz, can be used to sim- 



http://lambda.gsfc.nasa.gov/ 



* http://camb.info 
' http://class-code.net 
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Table 1. Summary of the main astrophysical components included in the current version of the model. See text for detailed descrip- 
tion of each component. 



Component 



Description 



Dipole WMAP dipole (Jarosik et al. 201 1), or user-defined. 

CMB Gaussian (constrained or unconstrained), non-Gaussianity, lensing component. 

Synchrotron Haslam et al. (1982) normalised and extrapolated with a power law (Miville-Deschenes et al. 2008). 

Free-free Dickinson et al. (2003) with normalisation of Miville-Deschenes et al. (2008). 

Spinning dust Draine & Lazarian (1999) with normalisation of Miville-Deschenes et al. (2008). 

Thermal dust Two-component model 7 of Finkbeiner et al. (1999). 

CO CO7=l-»0, 7=2->l, y=3->2 using Dame et al. (2001) and line ratios. 

SZ Cluster surveys, cluster counts, semi-analytic simulation or A'-body+hydrodynamical simulation. 

Radio sources Radio surveys extrapolated with power laws. 

Infrared sources IRAS survey modelled with modified blackbody emission laws. 

CIB Clustering properties of Negrello et al. (2004). 

UCH II regions IRAS catalogue with radio power law and modified blackbody function. 



ulate at random a population of clusters, and correspondingly, 
a model of SZ emission Pelabrouille, MeUn, & Bartlett|[2002| l. 
The same is true for point sources, on the basis of number counts 
as a function of flux density. 

Completely stochastic models of Galactic emission are more 
problematic, as the underlying statistical distribution of the emit- 
ting material is not known (if this concept makes sense at all). 
Nonetheless, a prescription to simulate realistic Galactic fore- 
ground emission is needed in any model of millimetre sky emis- 
sion. For such components, we adopt an approach which uses 
the predictive model based on existing observations, and com- 
plements it with stochastic corrections representative of current 
uncertainties in the data and the model. 

The generation of such constrained realisations, in which 
simulations match observations within observational eiTors, is 
implemented for most sky emission components in the present 
model (see, for example, paragraphs about CMB anisotropies in 
Sec. 3.2.3 galactic emission with small scales in Sec. 4.6 and 
constrained thermal SZ effect in Sec. |5.1.2] l. 



2.6. Use of the PSM in the Planck collaboration 

The PSM is the standard sky modelling tool used for simula- 
tions in the Planck collaboration. Maps and catalogues obtained 
from the PSM are either used directly in component separation 
separation challenges ( Leach et aT] |2008[ ), or used as inputs to 
the Level-S software for simulating Planck timelines (Reinecke 
et al. 2006 1. The PSM software package includes an interface 
with descriptions of the Planck HFI and LFI instruments for the 
generation of band-integrated sky maps, and the calculation of 
color-coiTection coefficients. 



avoids storing on local disks large data files which are rarely 
used (i.e., those that are used only for very specific values of in- 



2. 7. The model In practice 

In practice, the present sky model consists of a collection of data 
sets (maps and catalogues) and a software package, which can 
be used to produce maps of estimated or simulated sky emis- 
sion. The software presently available is written mostly in the 
IDL programming languagerland amounts to about 50,000 lines 
of code. Input data sets are stored in the form of ASCII and 
FITS files, for a total of about 500 GB of data. The software 
is designed to be usable on a modest computer (e.g., a stan- 
dard laptop); its data files are downloaded automatically from 
the associated data repository duiing the execution itself, which 



http://en.wikipedia.org/wiki/IDL_(programmingJanguage) 



put parameters). The HEALPix sky pixelisation scheme (Gorski 
|et al.| |200 5), in Galactic coordinates, is used in the current im- 
plementation. 

All the input parameters are passed using a single configu- 
ration file, which is read upon execution of the code. This con- 
figuration file is stored with the outputs of the simulation for 
traceability of the simulation. The software consistently uses one 
single master seed for the generation of all the necessary chains 
of random numbers, which guarantees the reproducibility of a 
set of simulations, and the independence of the various random 
numbers used for the generation of different components of sky 
emission and instrumental noise. 

The simulations can be generated at various resolutions, for 
various values of HEALPix pixel size, and for various values 
of maximum harmonic mode number {max- In principle, there 
is no limitation to the resolution (which can be set to 0), or to 
the HEALPix pixel size (which is limited only by the HEALPix 
software itself). The maximum CMB harmonic mode is set by 
the CAMB (or CLASS) software (reference CMB power spec- 
tra have been precomputed for {^-dx - 10,000). In practice, 
the model is constrained by observational data on scales of 6' 
to 1° for intensity, and few degrees for polarisation, depending 
on the component considered. The accuracy of extrapolations at 
smaller scales are expected to grow worse with increased res- 
olutions, and the model is expected to be useful down to the 
resolution of the Planck HFI (about 4' resolution). Test cases 
have been run up to HEALPix A?side=4096, and ^^ax - 6000. 
Technical information about the model, released software, data 
sets and documentation can be obtained from the PSM web site. 
PSM developers will make all possible efforts to maintain the 
code, document the consecutive versions, and give information 
about limitations and bugs on this website as well. 

The code runs in three consecutive steps: generation of a sky 
model; generation of sky emission maps at specific frequencies 
or in specified frequency bands, and generation of simulated data 
sets as observed by simulated instruments. 

2.7.L Step 1: the sky model 

The sky model is generated at a finite resolution that is set by the 
user and is common to all components. The sky model therefore 
implements maps of the sky smoothed with a Gaussian beam of 
full width at half maximum (FWHM) 0sky This ensures that all 
maps are near band-limited, and hence can be properly sampled 
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(with an appropriate HEALPix pixel size, which must be smaller 
than about one third of 0sky)- 

A particular sky model uses a common pixelisation scheme 
for all components (HEALPix, in Galactic coordinates, ring or- 
dering, and a given pixel size appropriate for sampling maps 
with no power above a specific harmonic mode /"max)- Polarised 
emission modelling is optional. 

The sky model consists of a set of maps of physical or 
phenomenological parameters describing the components: CMB 
temperature anisotropy; synchrotron amplitude and spectral in- 
dex maps; dust amplitude and temperature maps; catalogues of 
point sources (described by their coordinates and by their emis- 
sion laws for intensity and polarisation); catalogues of clusters 
(described by their coordinates, mass, redshift, temperature), etc. 
The sky model generated with the software is saved at the end of 
this first step; i.e., all the data sets used to model the sky are part 
of the output of a given run. It is then possible to apply steps two 
and three to the model at any time later on, without the need to 
re-run step one. 

2.7.2. Step 2: sky emission maps 

Sky emission maps result from the calculation of the observed 
emission at given frequencies, or in given frequency bands, on 
the basis of the sky model generated in the first step. The sky 
emission maps are what an ideal noiseless instrument would ob- 
serve with the specified frequency bands and with ideal Gaussian 
beams with FWHM equal to the resolution of the model sky flsky, 
and are sampled using the same HEALPix scheme as component 
maps. Note that this step is not fully independent of Step 3, as it 
is the list of instruments used to 'observe' the sky which sets the 
list of frequency bands used in Step 2. 

2.7.3. Step 3: observations 

Finally, sky emission maps generated in Step 2 can be 'observed' 
by simple model instruments. Each channel of the instrument 
is defined with a frequency band, a beam, polarisation proper- 
ties, an observing scheme, noise properties, and a format for 
the generated data sets. Although the generation of instrumen- 
tal noise excludes at the present slow drifts correlated between 
detectors, which generate signals that interfere with component 
separation (Delabrouille, Patanchon, & Audit 2002 Patanchon 



CMB dipole 



et al. 2008 1, it is possible to include the generation of such ef- 
fects in future developments of the code. 

The LEVEL-S software (Reine cke et al.[ |2006) can also be used 
as a means for sophisticated instrument simulation. 



2.8. Limitations 

The pre-launch Planck Sky Model is based on our understanding 
of the sky prior to Planck observations. It is intended to make 
available representative maps, to be used essentially in simula- 
tions. It should not be used for inferring the properties of the real 
sky. Limitations specific to some of the models of the individual 
components are specified at the end of each of the relevant sec- 
tions. 



3. Cosmic Microwave Background 

CMB anisotropies are modelled in the form of temperature fluc- 
tuations of a perfect blackbody. The assumed default CMB tem- 
perature is TcMB = 2.725 K ( |Fixsen et all[T996l|Fixsen|[2009l l. 




3354 uK_CME 



Fig. 1. CMB dipole prediction as generated by the code (here for 
a 1° 'sky resolution', in Galactic coordinates). 



Brightness fluctuations are modelled on the basis of a first order 
Taylor expansion of a blackbody spectrum. 



Icmb(0,<P) ^ Acmb(v)AT(0,,P), 



(1) 



where AT(6, (p) is a map of temperature anisotropies, and where 
the emission law A(v) is the derivative with respect to tempera- 
ture of the blackbody spectrum: 



iCMB 



(v) 



dT 



(2) 



T=Tr 



We ignore any effects, such as a non-zero chemical potential, 
which move the distribution of the CMB photons away from a 
perfect blackbody. 

The same emission law is used for CMB polarisation, which 
is represented with polarisation maps Q{9,^) and U{6,(p) or, 
alternatively, by the maps of scalar and pseudo-scalar modes, 
E(e, (f>) and Bje , 4>) (Za ldarriaga & Seljak|[T997l|Kamionko"wsE; 
[Kosowsky, & Stebbins|]T997] [ 

3.1. Tine CMB dipoie 

The CMB dipole is mostly due to the motion of the solar system 
with respect to the background. It is modelled on the basis of the 
measurement obtained with WMAP seven-year data (Jarosik et 
|aL]|201 1\ . The measured amplitude and direction are 

Adip = (3.355 ± 0.008) mK; 
Zdip = 263.99° ±0.14°; 

^7dip = 48.26° ±0.03°. (3) 

A dipole prediction uses by default these measured values. 
Alternatively, these values can be supplied by the user for set- 
ting up a different predicted dipole. 

A dipole simulation generates amplitude and direction as- 
suming a Gaussian distribution with default average values 
and standard deviations as given by the WMAP measurement. 
Alternatively, the simulation can use amplitude, direction and 
error as supplied in the input parameter list. 

Upon generation of a model sky, the amplitude of the CMB 
dipole map is modified to take into account the smoothing by 
the sky beam (even if for small beams the effect is negligible 
in practice). The default predicted CMB dipole generated by the 
PSM (here for a 0sky = 1° sky map resolution) is displayed in 
Fig.0 
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Fig. 2. CMB intensity and polarisation power spectra used by 



default (for Gaussian CMB simulations) - black: ; red: C 
(dashed parts are negative); green: C^^; blue: C^* (from scalar 
modes only). 
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3.2. Temperature and polarisation anisotropies 



On smaller scales, the CMB map is modelled as the outcome of 
a random process (the random generation of perturbations in the 
early universe). In the simplest model, fluctuations in the CMB 
are assumed to be Gaussian and stationary, so that its statistical 
distribution is entirely determined by its power spectrum C(. In 
addition to this generic model, the software implements a non- 
Gaussian CMB model, as well as predictions and constrained 
reaUsations matching WMAP observations. 



CMB anisotropies 

U Polarisation 



3.2.1. Gaussian CMB anisotropies 

Anisotropies on scales smaller than the dipole are imprinted on 
the CMB mostly at z ^ 1100, when primordial nuclei cap- 
tured free electrons to form neutral atoms. In the simplest model, 
these primary anisotropies are assumed Gaussian and stationary. 
Additional secondary perturbations arise due to electromagnetic 
and gravitational interactions, while photons propagate towards 
us. Reionisation of the universe at late times smoothes the small 
scale structure of the acoustic peaks, and generates 'reionisation 
bumps' in the polarisation power spectra. Gravitational interac- 
tion generates large-scale anisotropies via the integrated Sachs- 
Wolfe (ISW) effect, and modifies small scale anisotropies via 
lensing by large-scale structure (galaxies, and clusters of galax- 
ies). 

The statistics of the temperature and polarisation fluctuations 
of a Gaussian CMB are fully described by the multivariate tem- 
perature and polarisation power spectrum Cf of one tempera- 
ture map T and two polarisation maps E and B. For a standard 
cosmology (ignoring such effects as, e.g., cosmological birefrin- 
gence), the only non-vanishing terms in the power spectrum are 
Cj^, C[^, Cf^, Cf^, the other two (C[^ and C^^) vanishing for 
parity reasons. Fast and accurate software is available to com- 
pute these power spectra for a given cosmological model, and 
to randomly generate CMB temperature and polarisation maps 
for any given such spectrum: in the current version, interfaces 
to the CAMB ([Lewis, ChaUinor, & Lasenbyj [2000) l and CLASS 
( [Bias, Le sgourgues, & Tram||2011| l software are implemented 
to compute CMB power spectra for a given cosmological model. 
Alternatively, a default CMB power spectrum can be used, which 
matches current observational data. Fig. |2] displays the power 




Fig. 3. CMB temperature and polarisation anisotropy predictions 
as generated for a 1° beam, A^side = 256, in Galactic coordinates. 
This prediction is based on the CMB extracted from WMAP 5- 
year data using an ILC in needlet space. Note that the maps are 
not fully exempt from contamination by foregrounds and noise 
(see Del abrouille et al.] |2009| for a complete discussion of the 
component separation method used to obtain the CMB tempera- 
ture map). 



spectra used in the code for the current cosmological best fit, 
with no primordial 
of E modes only). 



with no primordial tensor modes (and hence C?^ from lensing 



3.2.2. CIVIB prediction 

CMB anisotropies have been observed by a number of experi- 
ments. In particular, CMB temperature maps have been obtained 
from publicly released WMAP data by a number of authors (e.g. 



Tegmark, de Oliveira-Costa & Hamilton j2003 ; Eriksen et al. 
|2004( iDelabrouille et al.| |2009| |Basak& Delabrouille) |20T2f . 
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Each of these maps can be used to predict the expected level of 
CMB anisotropy at any point of the sky. In the current version 
of the model, we use the Wiener-filtered Needlet ILC (NILC5) 
map obtained from WMAP 5-year data by |Delabrouille et aL] 
( |2009f ) as the standard CMB map at the resolution of the WMAP 
W-band channel. The harmonic coefficients x^^^ of that map are 
connected to the true CMB harmonic coefficients a^^^ by 



(4) 



where are the coefficients of the expansion of the sym- 
metrised beam of the WMAP W-band channel (as released with 
the WMAP 5-year data), we the coefficients of the harmonic 
Wiener filter applied to the NILC map, and «f„, is a residual ad- 
ditive term that accounts for all sources of error in the NILC5 
map (in particular instrumental noise and residual foreground 
emission). 

A predicted map at a different sky resolution (as specified in 
the list of input parameters of the code) is obtained from this one 
by multiplying the spherical harmonic coefficients by Bf/B*, 
where Bf are the coefficients of the expansion of the equivalent 
Gaussian beam corresponding to the target sky resolution 



-T _ T 



(5) 



Predicted CMB polarisation maps for Q and U are obtained as 
well, on the basis of the model TE correlation. The predicted 
E-modt harmonic coefficients are 



-E _ -T 



(6) 



where a^^^ are the predicted CMB temperature harmonic coeffi- 
cients computed from Equation [5] This polarisation prediction, 
however, is quite uncertain, since the TE correlation is weak. 
The predicted B-type polarisation vanishes in the current model 
implementation. Fig. [3] shows the predicted CMB temperature 
and polarisation at a target sky resolution of 1°. 

The predicted CMB maps are produced in such a way that 
the error term (difference between predicted map, and true CMB 
map at the target resolution) is minimised. They hence provide 
the best possible image of the true CMB sky given the present 
observations used in the model. With this criterion however, the 
power spectrum of the predicted CMB maps does not match the 
theoretical model. Indeed, given a noisy observation xg^ = s^m + 
riim, the Wiener filter as a function of I is, W[ - S [ jiS [ + N[), 
where S i is the power spectrum of the signal, and A^^ that of the 
noise. Then the power spectra of the predicted CMB temperature 
map is 

^TT l,_T |2\ ( r/2r^TT\ 



(7) 



i.e., the power spectrum is lower than expected, by a factor W( 
Similarly, the power spectrum of the predicted ii-mode polarisa- 
tion map is 

-EE ^ jCj'^) 
t £TT £EE 



[B]cr), 



(8) 



and the cross spectrum of T and E is 

-TE 



(9) 



3.2.3. Constrained Gaussian realisations 

Constrained CMB realisations, compatible with these observa- 
tions, are generated in the PSM on the basis of a theoretical 
model (power spectra of temperature and polarisation), and con- 
straints (the observed CMB temperature map described in 3.2.2 1. 

Recall that if {X, F)' is a two-dimensional random vector 
with mean (0, 0) and covariance matrix 



o-xY cr\ 



(10) 



then the conditional probability of X given Y is Gaussian with 
mean m and variance cp- given by 



and 



'Ay 



(11) 



(12) 



We simulate the constrained a^,,, by drawing random variables 
following the conditional distribution of a(,„ given some obser- 
vation "of,,,. 

The observation is performed at finite resolution, charac- 
terised by an effective beam B[. Each Of,,, of the CMB map is 
also measured with an error ef,„. The multipole coefficients of 
the observed map are 



(13) 



where af,„ and ef,„ are centred and independent Gaussian ran- 
dom variables, with variance Cc and A^^, respectively. Bf is the 
response of the beam in harmonic space (the beam is assumed to 
be symmetric). 

As cov(flft„,flf„,) = BfCf, var(flf,„) = C( and var(flf,„) = 
B^C; + N[ (we ignore correlations between the errors for dif- 
ferent (I, m), so that each mode can be generated independently 
of the others), we obtain, from the previous result, harmonic co- 
efficients of the constrained realisation that are Gaussian with 



BfCc 



B]Ce+Ne 



and variance 



CtNt 



BlCc + Ne' 



(14) 



(15) 



We note that as Nc — » oo or Bf — > 0, this law be- 
comes N(0,C[) (which is an unconstrained realisation of the 
CMB field), whereas if A^f - the conditional mean is 
Tiem/Bi. The mean value of the distribution of each harmonic 
coefficient matches the Wiener-filtered prediction described in 
Section llZ2l 



3.2.4. Non-Gaussian CIVIB anisotropies 

Non-Gaussian corrections to the statistics of CMB primary 
anisotropies are expected in the early Universe (inflationary) sce- 
narios. Software to generate non-Gaussian CMB maps has been 
dev eloped by|Smith & Zaldarriaga ( 2007 1, Liguori et al. ( 2007[ l 
and |Elsner & Wandelt| ( |2009l [ 



The first simulations of CMB temperature maps with primor- 
dial non-Gaussianity (NG) of the 'local type' have been based 
on the generation of the underlying primordial perturbation in 
Fourier space ( [Komatsu et al.[ |2003| l. A different method has 
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Linear CME map 



I 3S5 LiK_CMB 



Non-linear CMB correction (fp^ = 1) 




0.071 uK_CME 



Fig. 4. Linear part (top panel) and non-linear correction (bottom 
panel) for non-Gaussian CMB realisations used in the model. 
Note the different color scales, and that /nl ~ 1 (or a few) is 
a typical expected level. The cosmological parameters assumed 
for this simulation are from the fit of WMAP 7-year + BAO -H 
SN observations, as made available on the LAMBDA web site. 



then been proposed in Liguori et al. ( 2003 2007 1, where the 



authors work with filter functions to introduce the proper spa- 
tial correlations of the primordial potential. The current version 
of the model uses maps generated according to the method of 



Eisner & Wandelt (2009), a recent improvement of the method 
of [Liguori et al., ( ,2003| ) with better numerical efficiency. 

The non-Gaussian CMB model assumes a linear-plus- 
quadratic model for Bardeen's gauge-invariant curvature poten- 
tial oj^ namely 

<l)(x) = a)L(x) + /nl {<^L(xf - <<1)l(x)2)) , (16) 

where Ol is a Gaussian random field and the dimensionless non- 
linearity parameter /nl sets the level of NG in the primordial po- 
tential. The kind of NG described by this linear-plus-quadratic 
model is theoretically motivated by detailed second-order calcu- 
lations of the NG arising during — or soon after — inflation in 
the early Universe (Bartolo et al. 2004[ ); the assumption that the 
NG level is set by a constant parameter /nl is a fair approxima- 
tion, as long as |/nlI s> 1. 

Once realisations of the NG potential O are obtained, the 
harmonic coefficients for T and E are calculated by radial in- 



' Bardeen's gauge-invariant potential C) is defined in such a way that 
the temperature anisotropy reduces to AT/T = -0/3 in the pure Sachs- 
Wolfe limit. 



tegration of the linear and non-linear parts of the potential in- 
dependently, yielding maps of harmonic coefficients af„, for T 
and E and for both the linear and the non-linear part. The cor- 
responding input templates are used to generate simulated CMB 
skies for different values of /sjl. 

Because the computation of a non-Gaussian simulation at 
Planck resolution can take up to about 20 CPU hours, the code 
to make such simulations is not directly included in the pack- 
age. Instead, we use a number of pre-computed simulations with 
^max - 3500, and pick among those to generate a simulated non- 
gaussian CMB sky. 

Full-sky maps of one of the non-Gaussian simulations in- 
cluded in the model, for a 1 ° resolution simulation, are shown in 
Figure |4] where the top panel is the linear part and the bottom 
panel the non-linear correction (for fyi^ = 1). 

3.3. Lensed CMB 

One of the most important physical mechanisms that gener- 
ate secondary anisotropics is weak gravitational lensing of the 
CMB, which arises from the distortions induced in the geodesies 
of CMB photons by gradients in the gravitational matter poten- 
tial. As a result, the CMB temperature anisotropy that we ob- 
serve at a particular point on the sky in a direction n is coming 
from some other point on the last scattering surface in a dis- 
placed direction, such that. 



X(h) = X(n), 
h' — h + d{h). 



(17) 
(18) 



where X stands for any of the lensed CMB Stokes parameters /, 
Qor U and X is its unlensed counterpart. The directions h and 
h' are connected by the deflection field din). 

In the Born approximation, the deflection field is related 
to the Une-of-sight projection of the gravitational potential, 
'^(dAix)n,x) as d(h) - V<1)(«), such that. 



(D(n) 



1" 

Jo 



dx 



dhix^ -X) 



'V{dj,(x)h,x\ (19) 



where d^ix) is the comoving angular diameter distance corre- 
sponding to the comoving distance and Xs is the comoving 
distance to the last scattering surface. Although lensing depends 
upon all of the large-scale structure between us and the last scat- 
tering surface, most of the effect comes from structures of co- 
moving size of order few hundred Mpc at redshifts below 10. 
The typical deflection angle is around 2' to 3' but is correlated 
over several degrees. 

Weak gravitational lensing has several important effects on 
the CMB (see, e.g., Lewis & Challinor (2006i for a review). In 
addition to being an extremely robust prediction of the cosmo- 
logical concordance model, it is a probe of the structure forma- 
tion process in cosmology, overlapping with the onset of cos- 
mic acceleration, mainly in the linear or quasi-linear regime. 
The power spectrum of the deflection field is therefore useful for 
measuring parameters like the total neutrino mass (Lesgourgues 
l et al.j [2 006 ) or the dark energy equation of state (see, e.g. 
Acquaviva & Baccigalupi, 2006). 

Besides its intrinsic cosmological interest, the distortion of 
primary CMB anisotropics by lensing is a source of confu- 
sion for several scientific objectives of sensitive CMB experi- 
ments. Lensing modifies the CMB power spectrum, and gener- 
ates, through non-linear mode coupling, non-Gaussiarrity com- 
peting with the primordial one ( [Lesgourgues et al!||2005| l. It also 
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mixes primordial E and B modes, the main effect being the con- 
version of a fraction of the /i-modes into B-modes, which is a 
nuisance for the precise measurement of CMB B-modes from 
primordial gravitational waves. 

Weak lensing of CMB anisotropies by large-scale structure 
(LSS) has been measured on WMAP data by cross-correlation 
with the NRAO VLA Sky Su rvey (NVSS), used as a hi gh- 
redshift radio galaxy catalogue ( [Smith, Zahn, & Pore 2007 1, as 
well as in a combined analysis with SDDS and NVSS (Hirata 
et all |2008| l. Although marginally detectable in WMAP data, 
CMB lensing should be measured with high signal to noise 



by Planck without need to rely on an external data set (Hu & 
Okamoto 2002| l. However, in order to carry out this measure- 
ment in practice, the impact of instrumental (anisotropic beams, 
missing data, correlated noise) and astrophysical (Galaxy con- 
tamination, point sources, etc.) systematic effects on the CMB 
lensing estimators has to be studied with great care. This calls 
for the implementation, in the present model, of fast and accu- 
rate methods to simulate the lensing of CMB maps. 

In principle, such lensing should be compatible with the dis- 
tribution of matter between the last scattering surface and the ob- 
server This matter distribution is connected to the distribution of 
galaxies contributing to the anisotropies of the cosmic infrared 
background, to the radio source background, and to high-redshift 
galaxy clusters, which are simulated as well by the PSM. This 
connection between the different components of sky emission 
is ignored in the current PSM implementation, but will be the 
object of future improvements. 

3.3.1. Use of Non-equispaced FFT 

In order to compute the lensed CMB field on HEALPix grid 
points, we need to resample the unlensed CMB at arbitrary po- 
sitions on the sphere. The method implemented for doing so is 
based on the possibility to recast remapping on a sphere into 
remapping on a 2-d torus, and uses the non-equispaced fast 
Fourier transform (NFFT) to compute lensed CMB anisotropies 
at HEALPix grid points. The method is based on the previous 
work of Basak, Prunet, & Benabed (2009 ), which we briefly re- 
view here. 

The basic idea of the NFFT is to combine the standard fast 
Fourier transform and linear combinations of a window function 
that is well localised in both the spatial and frequency domains. 
Suppose we know a function / by means of evaluations f/^ 
in the frequency domain. According to the NFFT, the Fourier 
transform of that function evaluated at M non-equispaced grid 
points Xj in the spatial domain can be written as. 



k=-N/2 



■ m) X 



Inim k 
o-N 



fk 



^{2nkla-N) 



(20) 



where 0(^) is a window function in Fourier space, <I)(^) its coun- 
terpart in pixel space, and cr is the oversampling factor The in- 
dices j of the grid points Xj take the values j - 1 , 2, 3, . . . , M. 

The efficient evaluation of f{x) on irregularly spaced grid 
points requires a window function <1)(^) that is well localised 
both in space and frequency domain. The computational com- 
plexity is O {crN log N + K M), where K is the convolution 
length, cr is the oversampling factor, M is the number of real 
space samples, and is the number of Fourier modes. Among 
a number of window functions (Gaussian, B-spline, sinc-power. 



Kaiser-Bessel), the Kaiser-Bessel function is found to provide 
the most accurate results ( |Fourmontl |2003[ |Kunis & Potts| 
2008) , and is used in our simulations by default. 

3.3.2. Simulation procedure 

We summarise here the main steps of the simulation procedure 
of lensed CMB maps: 

1. Generate a realisation of harmonic coefficients of the CMB, 
both temperature and polarisation, and lensing potential 
from their respective theoretical angular power spectra ob- 
tained using CAMB; 

2. Transform the harmonic coefficients of the unlensed CMB 
fields and of the displacement field into their 2-d torus 
Fourier counterparts using the symmetry and recursion rela- 
tion of the Wigner rotation matrix (Varshal ovich, Moskalev] 
& Khersonski, 1988; Chalhnor & Lewis 200^, 

3. Sample the displacement field at HEALPix pixel centres us- 
ing the NFFT on the 2-d torus, apply this displacement field 
to obtain displaced positions on the sphere, and compute the 
additional rotation needed for the polarised fields using iden- 
tities of the spherical triangle (|Lewisl|2005l); 

4. Resample the temperature and polarisation fields at the dis- 
placed positions using the NFFT to obtain the simulated 
lensed CMB fields, sampled at the HEALPix pixel centres. 

In Fig. |5] we show one realisation of unlensed CMB tem- 
perature anisotropies, together with the lensing correction, the 
amplitude of the deflection field and the lensing potential map. 
We show a small portion of the same set of maps in Fig.|6]to il- 
lustrate the lensing effect more clearly. The correlation between 
the deflection field and the difference between the lensed and un- 
lensed CMB temperature anisotropies is clearly visible. We have 
found that, for unlensed CMB and deflection field maps with 
■^side - 1024 and /'^ax = 2048, 10"*^ accuracy is easily reached 
when setting the icr,K) parameters to (2,4). The current imple- 
mentation of the method to simulate lensed CMB with these 
parameters requires 15 minutes wall-clock computing time, us- 
ing 7.9 GB of RAM, on an ordinary PC configuration. Figure |7] 
shows that the average angular power spectra recovered from 
1000 realisations of lensed CMB maps (A^side - 1024) are con- 
sistent with the corresponding theoretical ones. This agreement 
between our numerical power spectra estimates and the CAMB 
predictions is both a validation of our numerical method, and of 
the CAMB estimates based on partially re-summed perturbative 
calculations. 

The accuracy of the lensing simulation is set by default in 
the current implementation, and is not a parameter of the soft- 
ware. If necessary it can be improved by increasing both the 
oversampling factor and the convolution length, at the expense 
of extra memory consumption and CPU time. Since the full pre- 
computation of the window function at each node in the spa- 
tial domain requires storage of large amount of real numbers 
{[2K+lf per pixel in 2-d), we have used a tensor product form 
for the multivariate window function (the default method within 
the NFFT library), which requires only unidimensional pre- 
computations. This method uses a medium amount of memory to 
store 2 [2K+1 ] real numbers per pixel, but at the cost of extra mul- 
tiplication operations to compute the multivariate window func- 
tion from the factors. In addition to that, pre-computation of win- 
dow function at harmonic domain requires storage of 2[2^'inax+2] 
real numbers. 

An increase in the convolution length not only increases the 
computational cost of the interpolation part of the NFFT, but also 
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Fig. 5. A simulated lensed CMB temperature anisotropy map (Top left), the amplitude of the difference of a simulated lensed and 
unlensed CMB temperature anisotropy map (Top right), the amplitude of the simulated deflection field map (bottom left) and the 
lensing potential map (bottom right) with HEALPix pixelisation parameter NsHe - 1024. These maps are obtained using the NFFT 
for the overs ampling factor cr -2 and convolution length K - A. 



increases the cost of the pre-computation of window function 
and memory requirement as one has to compute and store the 
window function at a larger number of grid points in the spatial 
domain before applying the NFFT. On the other hand, increas- 
ing the overs ampling factor only impacts the memory and CPU 
requirements of the (overs ampled) FFT part of the algorithm. 

The simulated displacement fields are computed from the 
gradients of Gaussian potential fields, neglecting non-linearities 
produced by the growth of structures and rotation induced by de- 
parture from the Born approximation. For CMB studies on scales 
comparable or larger than arc-minutes, this approximation is ex- 
cellent (see e.g. Bartelmann & Schneider) [2001) , and enables us 
to compare our lensed power spectrum estimates to the CAMB 
predictions, for which analytical estimates exist. However, our 
method is valid for any given displacement field; hence we could 
relax the Born approximation if needed, replacing it with ray- 



tracing in dark matter A^-body simulations ( Carbone et al. 2008 
2009[ l. Ray-tracing is affected by similar problems as the simu- 
lation of the lens effect on CMB maps, i.e., difficulties in accu- 
rately resampling a vector field on the sphere. Current state-of- 
the-art ray-tracing algorithms, such as described by [Teyssier et 
al. ( 2009) 1, could be made more accurate by using the technique 
described here. 



3.4. Limitations of the CMB models 

The various CMB models implemented in the PSM suffer from 
the following limitations: 



The prediction model is based on a CMB map extracted from 
WMAP data. That map is contaminated by noise and residu- 
als of foreground emission, particularly in the vicinity of the 
Galactic plane. It has been Wiener-filtered to minimize the 
total error. This however results in a filtered CMB map, the 
power spectrum of which hence is not that expected for the 
cosmological model considered. 

The CMB lensing and IS W as currently modeled are not con- 
nected to the distribution of large-scale structures. This will 
be developed in a future version. 

The non-Gaussian CMB is limited to ^n,ax - 3500, although 
the Gaussian signal can be extended to higher harmonic 
modes. 

Lensing of the CMB is currently implemented in the PSM 
pipeline only for Gaussian CMB models. It is possible, how- 
ever, to use the ilens code distributed with the PSM to lens 
Non-Gaussian CMB maps in a post-processing step. 



4. Diffuse Galactic emission 

Diffuse Galactic emission originates from the ISM of the Milky 
Way. The ISM is made up of cold atomic and molecular clouds, 
of a warm inter-cloud medium that is partly ionised, and of a hot 
ionised medium. The ISM also comprises magnetic fields and 
cosmic rays. Galactic emission from the ISM is usually sepa- 
rated into distinct components according to the physical emis- 
sion process: synchrotron radiation from relativistic free elec- 
trons spiralling in the Galactic magnetic field; free-free emission 
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Fig. 6. A portion of a simulated lensed CMB temperature anisotropy map (top left), the amplitude of the difference of a simulated 
lensed and unlensed CMB temperature anisotropy map (top right), the amplitude of the simulated deflection field map (bottom left) 
and the lensing potential map (bottom right) with HEALPix pixelisation parameter A^side = 1024. These maps are obtained using the 
NFFT for the oversampling factor cr -2 and convolution length K - A. The size of the maps displayed here is about 16° on a side. 



from the warm ionised medium, due to the interaction of free 
electrons with positively charged nuclei; thermal (vibrational) 
emission from interstellar dust grains heated by radiation from 
stars; spinning dust emission from the dipole moment of rotating 
dust grains; line emission from atoms and molecules. As the in- 
terstellar matter is concentrated in the Galactic plane, the inten- 
sity of these diffuse emissions decreases with Galactic latitude 
approximately according to a cosecant law (the optical depth of 
the emitting material scales proportionally to 1 / sin \b\). 

The Galactic diffuse emission carries important informa- 
tion on the cycle of interstellar matter in the Milky Way and 
its link with the star formation process. The precise mapping 
of the microwave sky by the COBE and WMAP experiments 
and the detailed analysis of the Galactic contamination to the 
CMB anisotropics have contribued significantly to our under- 
standing of the Galactic interstellar medium dWright et al.||1991| 
Bennett et al. 1992, Boulanger et al. , 1996; Finkbeiner, Davis, 
& Schle gel, ,1999, ,Banday et a H '2003; Bennett et al.^2003 



Finkbeiner| |2004[pavies et al.[ |2006; Mi ville-Deschenes et al 



2008). This has even led to the discovery of a new Galactic emis- 
sion component: the anomalous microwave emission (Kogutet 
aT] 1 1996a I. The Planck experiment will certainly expand our 



with its first series of papers (Planck Collaboration XIX.]|201 1 
Planck Collaboration XX 2011; Planck Collaboration XXI V 
|20riT|Planck Collaboration XXV||201 1) . 



In the context of the study of CMB anisotropics. Galactic 
emission is a nuisance. Its impact on the accuracy that can be 
reached on the cosmological parameters has been a concern 
for a long time and it is still paramount. For this reason, there 
have been significant efforts to model Galactic emission over the 
whole sky in the frequency range where the CMB fluctuations 
can be best observed. Many studies deal with specific Galactic 
components: sy nchrotron (|Giardino et aL 2002; Platani a et al. 
2003), free-free ([Dickinson, Davies, & Davis 2003 ; Finkbeiner 



2003 1, thermal dust ( Finkbeiner, Davis, & Schlegel 



more recently, polarized synchrotron and dust (jWaelkens et al. 





|2009t iFauvet et al.| [SoTT] jO'Dea et al.| |20I 2p. Some address 
the question more globally (Bouchet, Gispert &: Puget[ |1996[ 
iBouchet & Gispert[ ,1999; |Tegmark et al., ,2000 1, trying to esti- 
mate the foreground contamination as a function of frequency 
and scale. Because of the highly non-stationary statistical prop- 
erties of the Galactic diffuse emission, this exercise has limita- 
tions that can be partly overcome by using a tool that predicts 
maps of foreground emission for specific regions of the sky. 



1999|l and 



knowledge on the Galactic interstellar medium; it already has 
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Fig. 7. In all panels, the solid black line is the theoretical angu- 
lar power spectrum of the lensed CMB, and red filled circles are 
the average angular power spectrum recovered from 1000 reali- 
sations of lensed CMB maps at HEALPix A^sije = 1024. Lensed 
CMB maps are obtained using the NFFT for the overs ampling 
factor cr -2 and convolution length K - A. 



In that spirit de Oliveira-Costa et al. ( 2008| l developed a tool, 
based on a Principal Component Analysis decomposition of ra- 
dio survey maps, to produce maps of the total-power Galactic 
diffuse emission in the 10 MHz to 100 GHz range with an an- 
gular resolution of 1° to 5°. The PSM is similar but also models 
the sub-millimetre and far-infrared sky dominated by thermal 
dust emission, a dominant foreground in much of the Planck 
frequency range, and the polarised Galactic emission. The PSM 
also has the capability to produce simulations down to arcminute 



resolution. Also, contrary to the model of |de Ohveira-Costa et| 
al. ( 2008| l, the PSM is strongly linked to physical models of the 



Galactic emission processes. It is intended not only to simulate 
foreground emission for CMB studies but also to serve as a tool 
to gather relevant knowledge and models of the Galactic inter- 
stellar medium. 

For each Galactic emission component, the emission across 
frequencies is modeled using template maps (either observa- 
tions, or model maps that are constructed on the basis of existing 
data), which are interpolated in frequency according to specific 
emission laws. Each emission law depends on a set of physi- 
cal or empirical parameters, which can vary across the sky. Note 
that Galactic emission is actually implemented as a collection of 
different models using various input templates, and any of these 
can be easily called to produce different results. The PSM has 
been designed in such a way as to facilitate the integration of 
new models, in particular the ones coming out of the analysis of 
the Planck data themselves. At the time of writing, much has al- 
ready been learned about the Galactic emission with the Planck 
data. This knowledge is not in the PSM yet but it will be after 
the Planck data become public. 

In the following we describe the different models imple- 
mented for Galactic emission and in particular the combination 
of models which we believe, at the current time, is the most 
plausible and complete given the available data. The default 
Galactic emission model of the PSM is based on the combina- 
tion of the work of M iville-Deschenes et a l. (2008) who mod- 
eled the 23 GHz WMAP data by a sum of synchrotron, free-free 



and spinning dust emission, and the model of Finkbeiner, Davis, 
& Schlegel (1999) for thermal dust emission. In polarisation, 



the default PSM model includes 1) a synchrotron component 



based on the 7-year 23 GHz WMAP polarisation data ( Gold et 
al. 201 1} extrapolated in frequency with the same spectral in 



dex as for the total intensity, and 2) a dust component based on 



a new model presented in Sec. 4.7.2 This subset of models was 



chosen in order to reproduce the WMAP data, in intensity and in 
polarisation, while including an unpolarised spinning dust com- 
ponent. 

4.1. Synchrotron 

Galactic synchrotron emission arises from relativistic cosmic 
rays accelerated by the Galactic magnetic field (e.g.,'Rybicki &J 
[Lightman , 197 9 ). It is the dominant contaminant of the polarised 
CMB signal at low frequencies (below about 80 GHz). The in- 
tensity of the synchrotron emission depends on the cosmic ray 
density n^, and on the strength of the magnetic field perpendicu- 
lar to the line of sight. The spectral signature of the synchrotron 
emission is determined by the energy distribution of cosmic rays. 
In the simplest model, in the radio-mm frequency range, for elec- 
tron density following a power law of index p, {n^iE) oc E^''), 
the synchrotron frequency dependence is also a power law. 



^ sync 



(v) OC 



+2 



(21) 



or equivalently, in units of antenna (Rayleigh- Jeans) tempera- 
ture, 

rsync(v) oc (22) 

where /(v) and T(y) are related via the equation: 



/(v) 



2kv^ 



-T(v), 



(23) 



and where the spectral index, j3s - -ip + 3)/2, is equal to -3 for 
a typical value p - 3. 
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Synchrotron tomplote mop at 408 MHz 




Synchrotron tomplote mop at 408 MHz 




Fig. 8. Template map at 408 MHz used to model the intensity of 
synchrotron emission in the PSM. The colour scale of the bottom 
panel is histogram-equalised to increase the dynamic range. 



Synchrotron spectral index map 




Fig. 9. Synchrotron spectral index map used by default in the 
PSM. 



The synchrotron spectral index depends on cosmic ray prop- 
erties. It varies with the direction on the sky, and possibly, with 
the frequency of observation (see, e.g., [Strong, Moskalenko, & 



a general steepening of the spec trum towards high Galactic lati- 
tudes. In a more recent analysis, 'MiviUe-Deschenes et al. ( 2008| l 
find a lower scatter of the spectral index by taking into account 
the presence of additional emission in the radio domain from 
spinning dust grains. The spectral index map obtained in this 
way is consistent with ySj = -3.00 ± 0.06. Thus, as shown in this 
latter work, the inferred synchrotron spectral index map depends 
on assumptions made about other components. 

Due to energy loss of the cosmic rays, a break in the 
synchrotron spectrum is expected at frequency higher than ~ 
20 GHz. A detection of this effect has been reported by the 



WMAP team in the analysis of the WMAP first year data ( Bennett 
[et al., ,2003j . As mentioned in the analysis of the 7-year data, 
this apparent steepening could also be explained on the basis of 
a contribution of spinning dust to the total radio emission ( Gold 
[etal. 2011). 

Synchrotron emission is modelled on the basis of the tem- 



plate emission map observed at 408 MHz by Haslam et al. 
( 1982| l. The original map, reprojected on a HEALPix grid in 
Galactic coordinates at A^side = 512, and reprocessed to sup- 
press residual stripes and point sources, is obtained from the 
LAMBDA website. We further correct the map for an offset 
monopole of 8.33 K (Rayleigh-Jeans), which is subtracted to 



match the zero level of the template synchrotron of Giardino 
[et al.[ ( |2002l l. This offset includes the CMB (2.725 K), an 
isotropic background produced by unresolved extra-galactic ra- 
dio sources, the background monopole of synchrotron emission, 
and any zero level error in the data. 

This template synchrotron map is extrapolated in frequency 
using a spectral index map. Three options exist in the PSM for 
the spectral index. 



A constant value on the sky (set to -3.0 by default). 

The spectral index map of Giard ino et al.[ ( 2002| l obtained 

using a combination of 408, 1420 and 2326 MHz data. 



3. The spectral index map provided by model 4 of Miville- 
[Deschene s et al.[ ([2008) obtained using a combination of 
the 408 MHz and WMAP 23 GHz data. This model can be 
used with or without steepening at frequencies higher than 
23 GHz. 

This third model, without spectral curvature, is the default syn- 
chrotron model in the PSM. The template 408 MHz map used 
to model synchrotron emission is displayed in Fig. [8j and the 
spectral index map in Fig. 9] In the case of a simulation with 
resolution better than 1°, additional small scale fluctuations are 
added to the synchrotron intensity and spectral index maps, as 
described below in 14. 61 



4.2. Free-free 

Free-free emission is produced by electron-ion interactions in 
the ionised phase of the ISM. It is in general fainter than either 
the synchrotron or the thermal emission from dust, except in ac- 
tive star-forming regions in the Galactic plane. 

The free-free spectral index depends only slightly on the 
local value of the electronic temperature Tg and it is a slowly 



Ptuskin 2007 for a review of propagation and interaction pro 
cesses of cosmic rays in the Galaxy). 

Observations in the radio frequency range from 408 MHz to 



varying function of frequenc y (Bennett et al.[ 1992 Dickinson, 
Davies, & Davis[ |2003 ; Be nnett et al.[ [2003t . Between 10 and 
100 GHz, the spectral index varies only from 2.12 to 2.20 for 
4000 < Te < 10000 K. In the PSM, the spectral dependance 
of the free-free emission is based on the model described by 
10 GHz have shown spatial variations of the spectral index /3 Dickinson, Davies, & Davis (2003), assuming a constant elec- 



from -2.8 to -3.2 (Reich & Reich 


1988 Davies, Watson, &[ 


[Gutierrez[ 1 1 996| [Platania et al. 2003 1 


Bennett et al. 2003 1, with 



tronic temperature (Tg — 7000 K is the default value). The de- 
fault free-free emission law is plotted in Fig.[TO[ 
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Three options exists in the PSM to model the spatial varia- 
tions of the free-free emission. 



10 100 
Frequency [GHz] 



Fig. 10. Emission law of the free-free, spinning dust, and main 
CO molecular line emission (in units of spectral brightness ly, 
e.g., oc MJysr '). The spinning dust emission law is here nor- 
malised to unity at v = 23 GHz and the free-free emission law 
to /v = 2 at the same reference frequency. The integrated am- 
plitude of the first '^CO emission lines (transition (J=l-0)) is 
normalised to unity, and the other lines (transitions (J=2-l) and 
(J=3-2)) are displayed in proportion to their integrated intensity 
relative to the first one. 



Free-free map - 33 GHz 

Intenaity 




Free-free map - 33 GHz 

Intenaity 




Fig. 11. Template map at 23 GHz used to model the intensity 
of free-free emission in the PSM. As seen in the top panel, 
free-free emission is strongly concentrated in compact regions 
of the Galactic plane. The colour scale of the bottom panel is 
histogram-equalised to increase the dynamic range, and to show 
more extended, diffuse structures. 



1. The map of [Dickinson, Davies, & Davis ( |2003[ ) that uses a 
Ha emission template (a combination of WHAM ( Reynolds 
|etaLi[T998l ) and SHASSA dGaustad et all |200T] i data) cor- 
rected for extinction using the E{B — V) map of |SchlegeI^ 
[Finkbeiner, & Davis] ( |T998l ). 

2. The Maximum Entropy Method (MEM) free-free map of 
WMAP (Bennett et al.l 12003)). 

3. The free-free map obtained by Miville-Deschenes et al 



poos') that is a combination of the two previous maps. It uses 
the WMAP MEM decomposition of [Bennett et al.| ( |2003| ) in 
regions where the extinction is E{B - V) > 2 (or Ay > 6). 
In more diffuse regions, where Ha can be used more reli- 
ably as a proxy for free-free emission, it uses the estimate of 
[Dickinson, Davies, & Davis|p003) ), unless the WMAP MEM 
free-free is lower than the Ha estimate. 

These three maps have a resolution of 1°. The third one, dis- 
played in Fig. 11 is the default free-free model in the PSM. 



4.3. Molecular lines 

Molecular line emission is seen towards dense molecular clouds 
in our Galaxy and external galaxies. The strongest lines in 
the useful frequency range are those of '^CO at frequencies 
of 115.27GHz, 230.54GHz, and 345.80GHz (for (J=1^0), 
(J=2^1) and (J=3— >2) respectively). Other emission lines of in- 
terest, although fainter than those of '^CO, include: 

1. those of CO molecules involving other isotopes of carbon 
and oxygen (the most important being the '^CO emission at 
multiples of 110.20 GHz); 

2. a set of HCN lines around multiples of 88.63 GHz; 

3. the HCO+ lines at multiples of 89.19 GHz0 

Many other lines, such as those of CN, C2H and HNC, which 
have also been recently detected in high density molecular 
clouds in M82 (Naylor et al. 2010 1, are known to exist in the 
frequency range covered by the PSM. However, their contribu- 
tion to the observed broad-band emission in the frequency range 
covered by the PSM is small. 

The '^CO lines are of special interest, as they are the 
strongest ones located inside three of the Planck HFI frequency 
bands (the 100, 2 17, and 353 GHz channels; see Planck HFI[ 
[Core Team [201 1| ). They are thus important in the global model 
of sky emission. The '^CO emission is implemented in our 
model on the basis of the '"CO (J=l— >0) survey map of Dame, 



[Hartmann, & Thad deus (2001 ). The original map is not full sky, 
but contains most of the regions of strong emission (about 45% 
of the sky is covered, i.e., essentially all of the sky at Galactic 
latitudes \b\ < 30°). 

As surveys of other '^CO lines are missing, we choose to 
model the (J=2^1) and (J=3^2) lines by scaling the (J=1^0) 
map using the following line ratios in units of /TRjkms"' 
( ,Boulangerl[20T0| : 



(J=2- 
(J=3- 



>1)/(J=1- 
.2)/(J=2- 



>0) = 0.65; 
>1) = 0.50. 



The resulting relative integrated line intensities are displayed 
in Fig. 



10 



Emission from other molecular lines ('^CO lines of 
higher order or involving isotopes, and molecules other than CO) 
is neglected in the current model. 



http://www.splatalogue.net/ 
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4.4. Thermal dust emission 



Dust emission ratio 100/340 microns 



The thermal emission from heated dust grains is the dom- 
inant Galactic signal at frequencies above about 80 GHz 
Contributions from a wide range of grain sizes and composi- 
tions are required to explain the infrared spectrum of dust emis- 
sion from 3 yum t o 1 mm (Desert, Boulanger, & Puget 1990 Li 



& Draine 200 1[ Compiegne et al. 201 l[ l. At long wavelengths 
of interest for CMB observations (above 100 microns), the emis- 
sion from big dust grains at equilibrium with the interstellar ra- 
diation field should be the dominant source of the observed ra- 
diation. 

The thermal emission from interstellar dust is estimated on 
the basis of Model 7 of Finkbe iner, Davis, & Schlegel] ( |1999| l. 
This model, fitted to the FIRAS data (7° resolution), is based 
on the hypothesis that each line of sight can be represented by 
the sum of the emission of two dust populations, one cold and 
one warm. To limit the number of parameters, the model rests on 
three important assumptions: 

1. that the two dust populations are well mixed spatially and 
therefore both heated by the same radiation field; 

2. that the optical properties of the two populations are constant 
(i.e., fixed /3 and opacity e); 

3. that the ratio of their abundance (in mass) is constant. 

These assumptions imply that the ratio of absorbed (and 
emitted) power of the two populations is constant, a limita- 
tion that will be alleviated in the future using the 5' all-sky 
dust temperature derived from combining the IRAS and Planck 
data, such as computed injPlanck Collaboration XIX. (201111. At 



present, however, given these assumptions, Finkbeiner, Davis, 
|& Schlegel] ( |1999) ) showed that the spectral and spatial varia- 
tions observed are only determined by the local strength of the 
interstellar radiation field and the total dust column density. One 
can estimate the dust temperature of the two populations using 
only the ratio of the emission at two wavelengths. Following 



Finkbeiner, Davis, & Schlegel ( 1999 1, we use their map of the 



ratio of 100 jim to 240 jum emission based on DIRBE data, which 
has been processed to have a constant signal-to-noise ratio over 
the sky, resulting in a lower angular resolution (several degrees) 
at high Galactic latitude than in the plane (40')- This map is 
shown in Fig. 12 From the temperature of each dust popula- 
tion at a given position on the sky, the assumed ratio of emitted 
power, and the observed 100 fim brightness, one can estimate 
the dust column density A^, of the two populations at each sky 
position. The extrapolation at any frequency is then given by the 
sum of two modified blackbodies. 



(24) 



1=1 



where ByiJi) is the Planck function at temperature T,-, A^, is the 
column density of species /, and e,- v^' accounts for the normali- 
sation and frequency dependence of the emissivity. The 100 /zm 
map used is the one of [Schlegel, Finkbeiner, & Davis] ( |1998) l or 
a modified version of this map, displayed in Fig. 13 for which 
residual point sources, including a catalogue of ultra-compact 
Hn regions (see Section ]574] ), were removed. 

The main uncertainty in this representation of thermal dust 
emission comes from the model itself, rather than from noise in 
the observations used. We recall that the model is a fit to the 
very large scale dust emission (7°). Variations of dust proper- 
ties are observed everywhere in the ISM ( Planck Collaboration 




Fig. 12. Map of the ratio of 100/im to 240 jum emission, as ob- 
tained by ,Finkbeiner, Davis, & Schlegel ( 1999 1 and used in the 
present model. 

Tliermal dust emission 100— mioron template 




Fig. 13. Map of thermal dust emission 100 jum. Colours are 
saturated at lOOMJysr"', and the colour scale is histogram- 
equalised to increase the dynamic range and make both regions 
of low and high intensity visible. 



some models of interstellar dust predict significant variations of 
the modified blackbody spectrum for amorphous silicate grains 
(Boudet et al. 2005). The validity of the hypothesis of a constant 
ratio of absorbed power is therefore to be questioned. Under the 
hypothesis where the model used here is correct, the main uncer- 
tainty would come from the low angular resolution of the map 
of the ratio of 100 ;um to 240 //m emission which does not trace 
the small scale variations of the dust temperature. 



4.5. Spinning dust 

In the 10-100 GHz range several observations have pointed to 
excess emission, the so-called anomalous microwave emission 



(Kogut et al. 1996b 


Leitch et al. 1997| de Ohveira-Costa et 


al. ,1999 Watson et a 


. 2005|l. At frequencies above ~ 20 GHz 



]XXIV1 ]20Tll ]Planck Collaboration XXV] |2011| l. In addition. 



this difiiise emission has a spectrum falling with frequency, sim 
ilar to synchrotron, but has a spatial structure closer to the dis- 
tribution of thermal dust. One model proposed by Draine & 
[Lazarian] ( ]199 8) attributes this emission to small spinning dust 
particles. A detailed analysis of the WMAP data by Davies et 
al. ( 2006) 1 revealed the presence in selected areas of a com- 
ponent spatially correlated with thermal dust and with a spec- 
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tral signature compatible with the model of |Draine & Lazarian| 
1998| t. The spinning dust hypothesis is also compatible with 



free free ^ 1 dsQ r&s. 



[fgg fres with srr^ll ' 



3.les - 5 ercmin res. 



the analysis of the WMAP 23 GHz polarisation data (Miville- 
IDeschenes et al. , 2008 ), with the ARCADE 2 measurements be- 



tween 3 and 10 GHz ( .Kogut et al.||2011| l and with the fact that 
the microwave anomalous emission is well correlated with the 
12 fim emission divided by the radiation field strength ( Ysard,' 
Miville-Deschenes & Verstraete |201 0a), as predicted by models 
(Ysard & Verstraete 2010b| l. In addition, a recent joint analy- 
sis of data from the Planck mission, from WMAP, and from the 
COSMOSOMAS experiment shows unambiguous detection of 
a component in some specific regions, the spectrum of which 
peaks around 20 GHz and is compatible with spinning dust 
( [Planck Collaboration XX] [2011] ). 

Similarly to the free-free component, the spinning dust emis- 
sion is modelled using a spinning dust template map extrapo- 
lated in frequency using a single emission law. Three template 
maps are available. 



1. The E(B - V) map of [Sc"hlegel, Finkbeiner, & Dav"is|([T998]l. 

2. The spinning dust template based on model 2 of Miville- 
Deschenes et al. ([2008V This template was obtained us- 
ing the WMAP 23 GHz data and assuming a constant syn- 
chrotron spectral index of -3.0. 



3. The spinning dust template based on model 4 of Miville- 
[Deschenes et al.[ ( [2008] l. Same as the previous one but assum- 
ing a synchrotron spectral index constrained with the 23 GHz 
polarisation data. This model is the default one in the PSM. 

The emission law can be parameterised following the model of 
Draine & Lazarian| ( [l998[ l; the default adopted in the PSM, plot- 
ted in Fig. 10 is a mixture of 96% warm neutral medium and 4% 
reflection nebula. This specific parametrization of the [Draine &[ 
Lazarian[ ( |l998| l model has not been selected on any physical 
ground but because it fits the average spinning dust spectrum 
found in [Miville-Deschenes et al.[ ( [2008] ). 




1000 



2PO0 



30CW 



4.6. Adding small-scale fluctuations 

Small-scale fluctuations are automatically added to Galactic 
emission simulations when the resolution of the simulated map 
is greater than the resolution of the original template. The 
method used is described in [Miville^eschenes et al. (2007J . A 
Gaussian random field Gss having a power spectrum defined as 



^ [exp" 



• exp 



Icm J 



(25) 



is generated on the HEALPix sphere at the proper A^side- Here 
o"sinj and cr,g,„ are the resolutions (in radians) of the simulation 
and of the template to which small-scale fluctuations are added. 
The zero mean Gaussian field Gss is then normalised and multi- 
plied by the template map exponentiated to a power in order 
to generate the proper amount of small-scale fluctuations as a 
function of local intensity. This is motivated by the fact that the 
power of the Galactic emission fluctuations at a given scale gen- 
erally scales with the local average intensity at a given power 
between 2 and 3 (Miville-Desch enes et al.[|2007| l. The resulting 
template map with small-scale fluctuation added is then 



+ aGs,li 



(26) 



where a and /3 are estimated for each template in order to make 
sure that the power spectrum of the small-scale structure added 
is continuous with the large-scale part of the template. This ap- 
proach for generating small scales is similar to what has been 



Fig. 14. Top: a small patch of free-free emission before and af- 
ter adding random small scale fluctuations. Middle and bottom: 
power spectrum of free-free emission before and after adding 
small scales. 



done by Giardino et al. (2002 1, although in our case the weight- 
ing by the low-resolution map is slightly different. 

The addition of small scales is illustrated for free-free in 
Fig.[T4]and for spinning dust in Fig. [15] 

4.7. ISM Polarisation 

4.7.1. Synchrotron polarisation 

If the electron density follows a power law of index p, the syn- 
chrotron polarisation fraction is 



f,^3{p+l)/{3p + 7). 



(27) 



For p - 3, fs - 0.75, a polarisation fraction that varies slowly 
for small variations of p. Consequently, the intrinsic synchrotron 
polarisation fraction should be close to constant for a given vol- 
ume element in the Galaxy. However, once projected on the sky, 
geometric depolarisation arises due to variations along the line- 
of-sight of the angle between the magnetic field orientation and 
the line of sight direction. The magnetic field orientation de- 
pends on the spiral structure of the large-scale Galactic magnetic 
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Fig. 15. Top: a small patch of spinning dust emission before and 
after adding random small scale fluctuations. Middle and bot- 
tom: power spectrum of spinning dust emission before and after 
adding small scales. 



field and on local perturbation of the field due to turbulent mo- 
tions or shell expansions. The former can be estimated using a 
model of the spiral structure of the field ( Han et aL| 2006[ Sun 



etal. : 2008 i Sun & Reich 2010; Jans son et al.[|2009HJaffe et al. 
'2010 , 2011, Fauvet et al. , 2011), but the latte r can only be es- 
timated statistically ( [Miville-Deschenes et al.| [2008 ). Note that, 
since the synchrotron spectral index varies across the sky due to 
variations of the cosmic ray energy spectrum across the Galaxy, 
and since geometric depolarisation effects also depend on local 
effects in three dimensions, it is expected that the spectral index 
jSs is different for /, Q and U . We neglect this subtlety at this 
stage. 

Instead of relying directly on a model of the Galactic mag- 
netic field, our model of polarised synchrotron emission is built 
on the basis that synchrotron is the main low-frequency polarised 
emitter Hence we use the WMAP 23 GHz channel polarisation 
maps and U23, extrapolated in frequency using the same syn- 
chrotron spectral index as that of intensity. Other options may be 
implemented in the future, using either different models of the 
Galactic magnetic field, or interfacing the PSM with dedicated 




100 
Multipole, I 



Fig. 16. E and B power spectra of diffuse Galactic emission sim- 
ulated with the model at WMAP central frequencies (solid and 
dashed thick lines respectively). The P06 Galactic mask is used. 
The WMAP derived foreground levels from Gold et al. (2011) 
are also shown (thin lines). 



software for simulating Galactic synchrotron emission, such as 
the publicly available Hammurabi code (Waelkens et al. 2009). 

To extend the resolution of 223 and U23, small scales are 
added in a similar way as for the total intensity templates (see 
Eq.[26]l, but because Q and U are not strictly positive, the scaling 
of the Gaussian random field cannot be done with the template 
itself. Instead the small scale field is scaled with the unpolarized 



intensity field. For Q23 that is QL - Q23 + a Gj, li 



Fig. 



'23- 



17 displays synchrotron Q and U maps, as well as the 

The power 



polarisation fraction and angle as modelled here 
spectrum of the polarised sky emission model at WMAP fre- 
quencies, using the P06 mask, is displayed in Fig. 16 The figure 
illustrates the consistency of the simulated signal with WMAP 
results for polarisation ( Page et al.|[2007; Gold et al.||2011| l, at 
intermediate and high Galactic latitudes. It shows the Q for the 
sum of the polarised Galactic diffuse emissions evaluated in the 
WMAP K, Ka, Q and V bands, computed for simulation outputs 
and adopting the same WMAP P06 mask, excluding the bright- 
est Galactic emission. No beam convolution is applied. On large 
angular scales, a direct comparison with the spectra obtained by 
Gold et al. (2011 see also Fig. 17 in Page et al. 2007), shows 
reasonable consistency between the model and observed fore- 
ground levels. 

4.7.2. Dust polarisation 

Polarisation of starlight by dust grains indicates partial align- 
ment of elongated grains with the Galactic magnetic field (see 
|Lazarian [ |2007 , for a review of possible alignment mechanisms). 
Partial alignment of grains also results in polarisation of their 
emission in the far-infrared to millimetre frequency range. 

Very few observational data are currently available to model 
polarised dust emission, although things will drastically change 
soon with the observations of the Planck mission ( |Tauber et| 
al. 2010| . So far, dust polarisation measurements have mostly 



concentrated on specific regions of emission, with the excep- 



tion of the Archeops balloon-borne experiment (Benoit et al. 
12004 1, which has mapped the emission at 353 GHz over 17 



percent of the sky, showing a polarisation fraction around 4— 
5% in the Galactic plane, and up to 10-20% in some clouds. 
Recently, polarisation of dust emission at 100 and 150 GHz has 
been observed by BICEP over a smaller region of the Galactic 
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Fig. 17. Synchrotron Q and U maps at 23 GHz for a sky resolu- 
tion of 3°. At this frequency the synchrotron Q and U maps of 
the sky model are exactly the WMAP 7-year data. The two bot- 
tom panels show the polarisation fraction and the polarisation 
angle as implemented in our model. The average polarisation 
fraction is about 18%. 



pla ne (jMatsumura et al.| |2010^ and by QUaD ( |Culverhouse etj 



al. 2010 1. A large-scale estimate of the polarised dust power 



spectrum has been obtained with WMAP 7-year data by Gold 
jet al.j ( |201 l| l. The measured dust polarisation fraction is in 



rough agreement with what is expected from the polarisation of 
starlight (Fosalba et al.;'2002^. 

[braine & Fraisse^ ( ,2009) have shown that for particular mix- 
tures of dust grains, the intrinsic polarisation of the dust emis- 
sion could vary significantly with frequency in the 100-800 GHz 
range. As for synchrotron, geometrical depolarisation caused by 
integration along the line of sight also lowers the observed po- 
larisation fraction; the observed polarisation fraction of dust de- 
pends on its three-dimensional distribution, and on the geometry 
of the Galactic magnetic field. 

By making use of presently available data, we model po- 
larised thermal dust emission by extrapolating dust intensity, 
Ivip) (see Section 4.4 1, to polarisation intensity, assuming intrin- 
sic polarisation fractions, /j, for the two dust populations and a 
model for the dust geometric depolarisation, g^, and dust polari- 
sation angle, yj. The Stokes Q and U parameters for dust are 



Qv(p) = fd gdip) Ivip) cos(2yd(p)), 
Uvip) = fd gdip) Ivip) sin(27d(/?)). 



(28) 
(29) 



The value of is set to 0.15 for the two dust populations, con- 



sistent with maximum values observed by Archeops (Benoit et 



|ak 2004) and in good agi-eement with the WMAP 94 GHz mea- 
surement. It should be noted that does not correspond to 
the average observed polarisation fraction Pv/Iv, where Py = 
-y/f7|~+~2j. This is fdgdip), which has a mean of about 5% by 
default in the current model. Setting the same value of /d for 
the two dust populations implies that the polarisation fraction in 
the PSM does not change with frequency. The model has been 
designed to allow for different values of /d for the two dust pop- 
ulations, but in the absence of observational evidence, both dust 
populations have been given the same value for now. 

Because of the lack of observations, modelling polarised 
dust emission comes down to estimating maps of the geomet- 
ric depolarisation, g^, and of the polarisation angle, yd. One ap- 
proach is to use the large-scale structure of the Galactic magnetic 
field with a prescription for the turbulent depolarisation along 
the line-of-sight as O'Dea et al. (2012) did. This produces very 
smooth maps of ^d and yd as it does not reproduce the struc- 
ture projected on the sky due to the Galactic magnetic field as- 
sociated with the ISM turbulent dynamics. Such structure can 



be included by using a code like Hammurabi ( Waelkens et al. 
2009) that simulates the turbulent part of the magnetic field in 



three dimensions. This was done by F auvet et al. (2011) and is 
implemented as an option in the PSM but, even though the geo- 
metric depolarisation and polarisation angle maps produced that 
way have sensible morphology and statistical properties, their 
detailed structure does not match the observed sky. 

In order to simulate dust polarisation maps that stay as close 
as possible to the data, and to be coherent with the polarised 
synchrotron model, we make use of the 23 GHz WMAP polar- 
isation data to estimate the geometric depolarisation and polar- 
isation angle maps for dust. Here we explicitely make the as- 
sumption that, on average, synchrotron and thermal dust polar- 
ized emission are influenced by the same magnetic field and 
that their maps of g and y should be correlated. Therefore, in 
our model ^d and yd are based on equivalent maps estimated 
for the synchrotron emission using the WMAP 23 GHz and the 
Haslam 408 MHz data, the synchrotron polarisation angle 



ys = ^ tan ' (-t/as, Q23) , 



(30) 
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and geometric depolarisation factor 



g synchrotron 



1/2 



/,/408 (23/0.408) 



-3 ' 



(31) 



where all templates are smoothed to 3°. In order to imprint a 
specific signature of dust on and js, we have modelled what 
would be the large-scale structure of those two quantities for dust 
and synchrotron. In order to do that, we have built maps of /, Q 
and U for dust and synchrotron given a model of the Galactic 
magnetic field and a simplistic model of the distribution of dust 
and electrons in the Milky Way. The maps of g and y are then 
obtained by inverting Equations ( [28] l and (|29]l. 

Assuming only one dust population of volume density /id, 
the total intensity of dust is 

/d(v) = J«dedi^^Bv(rd)&, (32) 

where the integral is along the line of sight z- For Q and U the 
integrals include terms related to the angle between the magnetic 
field direction and the line of sight; 



and 



where 



cos 20 sin a t/z, 
Ui{v) - fd Jn^edV^^ByiTd) sin 20 sin a t/z, 



cos 2(p - 
sin 2<p = 



-2B,By 



sin a = ^1 -BlIB^, 



and 



(33) 

(34) 

(35) 
(36) 
(37) 

(38) 



with the X - y plane corresponding to the plane of the sky. For 
synchrotron, the unpolarised intensity / is given by 



/s(v) 



,(l+rt/2 



dz. 



(39) 



where Ps is the synchrotron spectral index defined in 
Equation (22 1. As in the dust case, the Stokes parameters of po- 
larised synchrotron emission are 



(40) 



and 



Qs(y) - fs J'^sV^'iieB^l^''^^^ cos 2(f> sin a dz, 

U,(v) = /, j^ey-n^B^'^'''^^ sin 20 sin a dz. (41) 

The polarisation fraction /, is related to the cosmic ray energy 
distribution according to Equation ( [27| . 

Equations from (32i to (41 1 show the difficulties of mod- 



elling polarised emission. This requires a three-dimensional 
model of all dust properties (rid, fd, /3d and Td), of the energetic 
electrons, and of the magnetic field geometry. To simplify this 
task we assume that £d, Pd and Td are constant along a given line 
of sight. Then, following |Page et al.| ( |2007| , |MiviUe-Deschenes| 



Y synchrotron (rad) 



Y dust (rad) 



0.46 



0.23 





0.00 



1.57 

-0.00 

-1.57 
1.57 

0.00 

-1.57 



Fig. 18. Maps of the large-scale geometrical depolarisation g' 
and polarisation angle y' for the synchrotron and dust, based on 
a model of the Galactic magnetic field and density distributions 
of the energetic electrons and dust grains. These maps are used 
to correct the dust depolarisation and polarisation angle maps 
deduced from the 23 GHz polarisation data for the difference 
in the large-scale structure between synchrotron and dust (see 
Equations [43] and|44l). 




et al. (2008 1 and Fauvet et al. (201 1 1, we model the density dis 



tribution of dust and electrons with a galacto-centric radius and 
height function 



n — HQ exp{-r / hr) sech (z/ hj, 



(42) 



with scaling parameters for the galacto-centric radius and height 
dependence equal to hr = 3 kpc and = 0.1 kpc for dust, and 
hr - 5 kpc and /i, = 1 kpc for electrons. 
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Finally we use the model of the magnetic field described in 
[Miville-Deschenes et al. ( 2008 1, which has a regular component 
following a bi-symmetrical spiral with a pitch angle of -8.5° 
and a turbulent component with an amplitude of 1.7 yuG (i.e., 
0.57 times the local value of the magnetic field Bq = 3 ^G). The 
turbulent part of the field used here cannot reproduce the multi- 
scale properties of the interstellar emission on the plane of the 
sky; it is included to take into account the fact that longer lines of 
sight in the Galactic plane have a stronger depolarisation effect, 
due to a larger number of turbulent fluctuations of the magnetic 
field than at high Galactic latitude. The turbulent part of each 
line of sight is simulated independently of its neighbours, so this 
method produces no turbulent structures on the plane of the sky. 

The maps of g and y obtained using this model reflect po- 
larisation effects on large scales (> 20°). These maps, shown in 
Fig. 



18 are similar to the ones used in the model of O'Dea et al.l 



Q (mKCMB) 



0.03 



2012| l. As mentioned previously, local variations of the dust and 



electron densities and of the magnetic field properties caused by 
the complex physics of the interstellar medium will produce spe- 
cific structures in the real maps that cannot be reproduced in this 
simple model. The objective here is to capture large-scale differ- 
ences between the synchrotron and dust polarised emission due 
to the fact that 1) dust grains and energetic electrons have dif- 
ferent scale-heights and lengths in the Galaxy; and that 2) their 
respective polarised emissions do not depend in the same way 
on the magnetic field (unlike the dust, the synchrotron emissiv- 
ity depends only on - see Eq. 39 1. Fig. 18 shows indeed that 
the geometrical depolarisation and polarisation angle maps of 
these two processes are slightly different. 

To produce the dust geometrical depolarisation {g^) and po- 
larisation angle (yj) maps we combine this large-scale model 
(that provides the structure at scales larger than 20°) with the 
synchrotron and js estimated with the 23 GHz and 408 MHz 



templates (see equations 30 and 31 1 to add structures between 



3° and 20°. The dust geometrical depolarisation and polarisation 
angle maps used in the model are then defined as 



and 



g'd + igs - g's). 



(43) 



(44) 



where g'. and y' are the maps built with the large-scale model and 



shown in Fig. 18 The addition of fluctuations at scales smaller 
than 3° is performed in a similar way as for the synchrotron po- 
larisation: they are added in Q and U directly with a scaling 
following the intensity /. 

The dust polarisation maps estimated with our model at 
200 GHz (in the frequency range of interest for high resolution 
CMB observations) are shown in Fig. 19 As expected, the po- 



larisation fraction and polarisation angle maps are very similar 
to the ones for the synchrotron shown in Fig. [17] On the other 
hand, the maps of Q and U are significantly different from those 
for the synchrotron, because of the distinct spatial structure of 
the thermal dust and synchrotron unpolarised emission. We ex- 
pect the modelling of polarised dust to be significantly improved 
in the future, with the use of the Planck polarised data between 
100 and 353 GHz. 



4.7.3. Polarisation of otiier ISM emission 



Free-free emission is intrinsically unpolarised. At second or- 
der, one expects some polarisation towards ionised regions by 
a mechanism similar to that giving rise to CMB polarisation: if 
incident light, scattered by free electrons in the ionised gas has a 
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Fig. 19. Model of the dust polarisation at 200 GHz and at 3° reso- 
lution. The polarisation fraction and polarisation angle maps are 
slightly modified versions of the same maps obtained for syn- 
chrotron at 23 GHz using the WMAP and 408 MHz data. 



quadrupole moment perpendicular to the line of sight, then scat- 
tered light is preferably polarised perpendicular to the direction 
of incidence of the strongest incoming radiation, by reason of 
the polarisation dependence of the Thomson cross section. This 
is neglected in the current version of the model. The Thomson 
scattered free-free light is expected to be small (a few percent or 
likely less, [Page et al.| |2007| [Gold et al.| |20n| [Macellari et al] 
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|201 with a maximum of about 10% on the edge of bright H ii 
regions (e.g., Rybicki & Lightmanl|1979| l. 

Spinning dust emission is observed to be at most weakly po- 
larised, of order a few percent (jBattistell i et al.j 2006 Lopez- 
ICaraballo etaTHZOTTllDickinson, Pee l, & Vidalf |20l7] l. On the 



theoretical side, Lazarian & Draine ( 2000) predict very low (but 
possibly measurable) levels of polarisation for spinning dust 
emission. Even if anomalous dust emission is due to magnetic 
dust (Draine & Lazarian 1999| l, the polarisation is expected to 
be low (~0.5%) if the iron inclusions are random ( [Lazarian &| 
Draine||20001l. In the current software implementation, given the 



absence of observational evidence or theoretical predictions of 
significant spinning dust polarisation, its emission is assumed 
unpolarised. 



4.7.4. Limitations of tine models of galactic emission 

There is no theoretical statistical model for ISM emission. 
Models for the emission of all the galactic components are 
phenomenological descriptions based on observed emission. As 
such, the models are at least as uncertain as sky observations. 

In particular, free-free emission, which does not dominate 
at any frequency, is poorly measured so far Maps of free-free 
emission used in the PSM are derived from either or both of an 
Ha emission composite map that combines two different sur- 
veys, and a map based on the WMAP MEM decomposition. The 
former, using Ha as a proxy for free-free emission, is subject to 
various uncertainties (dust absorption, scattering of Ha by dust 
grains, variations in electron temperature). The latter strongly 
depends on the model assumed in the WMAP analysis, and does 
not seem to match well the radio recombination line measure- 
Alves et al. ( 2010| l (albeit over a limited sky region). 



ments m 

Hence none of the two basic free-free templates used in the PSM 
is very reliable, and the modelled PSM free-free emission should 
be considered at best as indicative. 

Thermal dust emission is modeled following the study of 
Finkbeiner, Davis, & Schlegel ( 1999 1 based on an analysis of 
the FIRAS data. These authors showed that the data can be de- 
composed by the sum of the emission from two dust populations. 
First, it is important to note that, even though it is thought that 
interstellar dust is made of carbonaceous and silicate grains, the 
dust spectral indices and equilibrium temperatures used in this 
model are not compatible with the current knowledge of ther- 
mal emission from big grains (Li & Draine 200 1[ |Compiegnel 
|et al.| |201 1[ ). Combined with the limited signal-to-noise ratio 
of the FIRAS low frequency data used to develop this model, 
the spectral shape of the dust emission at frequencies lower than 
300 GHz is not expected to be very accurate. Second, the tem- 
plate of the dust temperature used in this model is based on 
a version of the 1 00/240 //m ratio map from the DIRBE data, 
smoothed at high Galactic latitude to increase the signal-to-noise 
ratio. Thus the effective resolution of the dust temperature tem- 
plate is varying across the sky, from 30 arcmin in the plane to 
only several degrees at high Galactic latitude. That has to be 
kept in mind when investigating the variations of the dust emis- 
sion properties at high Galactic latitude. Similarly, the template 
and emission laws of spinning dust emission are representative, 
but not exact. The models for thermal and spinning dust emis- 
sions will be adjusted in future releases of the PSM, after the 
analysis of the Planck observations has been completed. 

As discussed by |Sun et al.| ( |2008| , the synchrotron emission 
is absorbed at low frequency by foreground thermal electrons 
that produce free-free emission. This absorption is less impor- 
tant at 23 GHz than at 408 MHz, which produces an artificial 



flattening of the synchrotron spectral index. Therefore, [Sun et 
(2008 1 showed that it is likely that in the inner portion of the 
Galactic plane, where the free-free absorption can be significant, 
the 408 MHz map underestimates the synchrotron emission if 
extrapolated at higher frequencies with a constant spectral in- 
dex. At this point, given the significant uncertainty on the free- 
free emission in the Galactic plane and on the filling factor of the 
diffuse ionised gas (the opacity is proportional to the square of 
the electron density), it is difficult to quantify this effect on the 
408 MHz map. In addition, the fiducial synchrotron model of the 
PSM is not based on an extrapolation of the 408 MHz map; this 
synchrotron template has been derived from the WMAP 23 GHz 
polarisation data, corrected for a model of the synchrotron depo- 
larisation jMiville-Deschenes et al.|p008| l. In this case the uncer- 
tainty on the synchrotron map is coming from the model of the 
depolarisation itself. 



5. Emission from compact sources 

Numerous galactic and extragalactic compact sources emit radi- 
ation in the frequency range modelled in the PSM. Extragalactic 
emission arises from a large number of radio and far-infrared 
galaxies, as well as clusters of galaxies. Compact regions in 
our own galaxy include molecular clouds, supemovae remnants, 
compact Hn regions, as well as galactic radio sources. 

The thermal and kinetic SZ effects, due to the inverse 
Compton scattering of CMB photons off free electrons in ionised 
media, are of special interest for cosmology. These effects oc- 
cur, in particular, towards clusters of galaxies, which are known 
to contain hot (few keV) electron gas. SZ emission is modelled 
here both on the basis of known existing clusters, observed in 
the optical and/or in X-rays, and on the basis of cluster number 
counts predicted in the framework of a particular cosmological 
model. 

Distant radio and far-infrared sources are not resolved by 
CMB instruments, and are represented in our model as a pop- 
ulation of point sources, each of which is described by a number 
of parameters (type of the source, position on the sky, flux and 
polarisation as a function of frequency). Although unresolved, 
the brightest objects are detected individually in observations, 
and represented in the form of a catalogue (rather than a map). 

In addition to those brightest sources, a diffuse source back- 
ground arises from the integrated emission of a large number of 
objects that are too faint to be detected individually, but which 
contribute to sky background inhomogeneities. Background 
sources are represented using maps of resulting brightness fluc- 
tuations. 

Finally, our Galaxy itself contains a number of compact re- 
gions of emission that can be modelled as point sources. Some 
of those are modelled as part of the catalogues of radio and far- 
infrared sources (without distinction between galactic and ex- 
tragalactic objects). A specific population, identified as ultra- 
compact galactic Hii regions, is singled out and treated sepa- 
rately. 



5.1. Sunyaev-Zeldovich effects 

Hot electron gas, as found in the largest scale structures in the 
universe, interacts with incoming photons by inverse Compton 
scattering. In particular, background photons are scattered by hot 
gas in clusters of galaxies, giving rise to secondary anisotropics 
of the CMB by shifting the energy of the interacting photons. 
The general properties of the effect have been described first 
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by|Sunyaev & Zeldovich] (|T972li (see also |Birkinshaw| [T999| 
Carlstrom, Holder, & Reese 2002 for recent comprehensive re- 
views). 

Classically, a distinction is made between two SZ effects. 
The thermal SZ effect (tSZ) corresponds to the interaction of 
the CMB with a hot, thermalised electron gas, while kinetic SZ 
(kSZ) corresponds to CMB interaction with electrons having a 
net ensemble peculiar velocity along the line of sight. 

In the thermal case, the electron population is characterised 
by a Fermi distribution at a temperature Tg much larger than 
the CMB temperature (in typical clusters of galaxies, is of 
order few keV). For non-relativistic electrons, the change in sky 
brightness is; 

6Iy = yf(v)By(TcMB), (45) 

where Bv(7'cmb) is the CMB blackbody spectrum, /(v) is a uni- 
versal function of frequency that does not depend on the physical 
parameters of the electron population (and hence on the galaxy 
cluster they are associated with), and y, the Compton parameter, 
is proportional to the integral along the line of sight of the elec- 
tronic density rte multiplied by the temperature of the electron 
gas: 

/kT 
%rn^a-jdl. (46) 

For a standard cluster detectable at millimetre wavelengths from 
its thermal SZ effect with current instruments, - 5 keV, the 
optical depth of the gas, r = ^ n^a-jdl, is of order 10"^, y - 10""* 
and the Compton parameter integrated over the cluster angular 
extent, Fsz, is of order 10""^ arcmin^ for a typical cluster angular 
size 1'. 

In the kinetic case, the interaction of CMB photons with a 
population of moving electrons results in a shift of the photon 
distribution seen by an observer on Earth: 



6L 



dBy{T) 



dT 



(47) 



where /3r - vv/c is the dimensionless cluster velocity along the 
line of sight. For typical radial velocities of 300 km/s (J3r - 
IQ-^), and for r ^ 10"^ the net kinetic SZ effect is AT/T ^ 10"^ 
At present, only limited SZ observations are available. Those 
observations aim both at studying previously known clusters 



(|Basu et al.| |2010t |Zwart et al.| |2011| |Hincks et al. 2010 


and at discovering new objects blindly (jVanderlinde et al. 


2010llWilliamson et al.| 


20 11 , .Planck Collaboration VIII| |20 1 1 


Marriage et al.| 201 la| 


Reichardt et al.||2012i. They have not 



only confirmed the existence of the effect but they have also 
made possible a variety of studies which allow predictions of 
the general properties of the SZ sky, such as the shape and nor- 
malisation of its angular power spectrum ( jShaw et al.| |2009| l, 
and the connection between SZ observables and cluster masses 
or the X-ray luminosity from direct free-free emission of the hot 
ionised gas (jMelin et al!] [2011] [Planck Collaboration X| [20TT] 
IPlanck Collaboration XI| [20Tlt . 

Several options are available to simulate SZ emission with 
the present sky model. A prediction model generates ther- 
mal SZ effect as expected from catalogues of known clus- 
ters, observed in X rays with ROSAT or in the visible with 
SDSS (Section |5.1.1| l. A simulation of both thermal and ki- 
netic SZ emission is possible on the basis of cluster number 



counts from cluster mass functions (Section 5.1.2i, or using A^- 
body and hydrodynamical simulations of large-scale structures 
(Section 5.1.3i, or making use of hydrodynamical simulations 
at low redshift (z < 0.25), complemented by cluster counts at 
higher redshift (Section [5.1.4| i. 




Fig. 20. Density of MaxBCG clusters. The clusters are mainly 
concentrated in the northern Galactic hemisphere. 

5.1.1. SZ emission prediction 

For any given cluster of known mass and redshift, scaling re- 
lations exist that allow one to predict the level of thermal SZ 
emission originating from that particular cluster. The prediction 
of the thermal SZ emission of the sky in our model is based on 
the observations of 1743 galaxy clusters with the ROSAT X-ray 
satellite (NORAS ( Bohringe rltaL] |2000| , REFLEX (Bohring^] 



et al. 2004) , serendipitous surveys), and of 13,823 optically se- 



lected clusters extracted from the SDSS galaxy survey. 215 clus- 
ters are common to these two independent catalogues. 

The ROSAT catalogue used in the model is the Meta- 
Catalogue of X-ray detected Clusters of galaxies (MCXC) by 
|Piffaretti et al.| ( p011| l. It covers the full sky, although with un- 
even depth. These clusters are located at redshifts ranging from 
0.003 to 1.26, the median redshift of the sample being 0.14. For 
each cluster observed in X-rays, the predicted SZ flux Fjoo is 



given by the combination of the Lx-M relation of Pratt et al. 
( |2009l l and the M-Y relation of |Arnaud et"aL] ( [20T0l l. The cluster 
model adapted here and its normalisation on the basis of cluster 
observations is detailed in Planc k Collaboration Xl l'^Oll). We 
simply cite here the two key scaling laws used in the model. The 
first one connects the X-ray luminosity and the mass: 



E{z) 



-7/3 



L5OO 



1 044 erg s- 



M50Q 



3 X 1014 M, 



(48) 



where E{z) is the Hubble parameter normalised to its present 
value, L500 is the X-ray luminosity emitted, in the [0.1-2.4] keV 
band, from within the radius /?5oo at which the mean mass den- 
sity is 500 times the critical density of the universe at the cluster 
redshift, M^qq is the cluster mass within Rsoo, and the numerical 
values of the constants are log(Cz,M) - 0.274 and o^m = 1 -64. 
The second one connects the SZ flux and the mass: 



Fsoo = 1.383 x 10^^/(1) 



M500 



X E(z) 



2/3 



500 Mpc 



arcmin 



(49) 



where Fsoo is the SZ flux in Rsoo, 1(1) - 0.6145 is a numerical 
factor arising from the volume integral of the pressure profile, 
Z)a(z) is the angular distance and omyx - 0.561. 

MaxBCG clusters (Koester et al. 2007) are detected mostly 
in the northern Galactic hemisphere (see Figure 20 1. These clus- 
ters are located at redshifts ranging from 0.1 to 0.3, the median 
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whose mean is given by a mass function: Press-Schechte r (|Press| 
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Fig. 21. Top panel: Cluster number counts for different mass 
functions. Bottom panel: Cluster number counts for different val- 
ues of erg, using the Tinker mass function. 



redshift of the sample being 0.23. We compute the SZ flux F500 



using the Fsoo - A^2oo relation derived by the Planck Collaboration 

pail ( |2on] i, 



7500 - F20 E ' (z) arcmm , 

^"l 20 / \500Mpc/ 



(50) 



where A^2{)0 is the cluster richness, F20 - 7.4 lO^^arcmin and 
Q'20 = 2.03. 

The SZ prediction implemented in the current version of our 
model does not include kinetic or polarised SZ effects. 

5.1.2. SZ emission from cluster mass functions 

An approach for simulating the SZ effect from galaxy clusters 
consists in assuming a cluster mass function N{M,z)- Following 
the simulation method described in [Delabrouille, M elin, &| 
[Bartlett (2002), we start from a set of cosmological parame- 
ters (Hubble parameter h; density parameters for matter Q„, 
vacuum energy Qa, and baryons Qb; amplitude erg of fluctu- 
ations on 8/;"' Mpc scales). This set of cosmological parame- 
ters matches what is used for simulating CMB anisotropics (and 
hence is defined uniquely for both). The cluster distribution in 
the mass-redshift plane {M,z) is then drawn from a Poisson law 



|& Schechterl [T974|), S heth-Tormen dSheth & TOTmenl |1999| l. 



Jenkins Penkins et al.| 2001), Evrard ( Evrard et al. 2002), or 



Tinker ( [Tinker et al. 2008 ). Cluster Galactic coordinates (/, b) 
are then drawn at random on the sphere, with a uniform proba- 
bility distribution function (which here neglects correlations in 
the large-scale structures traced by galaxy clusters). Baryons are 
pasted on clusters assuming an electron density profile (either a 



Navarro, Frenk, & White (1996i profile, modified according to 
Nagai et al.| ( ,2007i ) and |Arnaud et alT|(|2010), or a j6-model ) and 



using a mass-temperature relation (Pierpaoh et al. 2003[ ), the 



normalisation of which is a free parameter, or a Ysz-M relation 



from Chandra orXMM observations ( .Nagai et aL]|2007[|Arnaud 
|etaL||20T0l ). 

The three-dimensional velocities of the clusters are drawn 
from a Gaussian distribution with zero mean and standard de- 
viation given by the power spectrum of density fluctuations 
( |Peebles||1993) ). This, again, neglects correlations between clus- 
ter motions, i.e., largest-scale bulk flows. The kinetic SZ map is 
obtained using this velocity field and the electron density of the 
clusters. 

This semi-analytic method of simulating SZ clusters requires 
relatively small amount of CPU time. It is hence fast enough to 
allow for a new generation of the cluster catalogue with each run 
of the code, and thus makes it possible to produce many differ- 
ent sky maps with different input parameters as, for example, a 
varying value of erg. The gas properties and scaling relations can 
be easily tuned to the observed ones (and match those measured 
recently from WMAP and Planck data). It is possible to include 
the contribution from clusters in an arbitrarily large mass range, 
including groups of galaxies. Simulations can be repeated many 
times with different seeds, which permits us to characterise the 



selection function of a data processing chain ( Melin, Bartlett, & 
Delabrouille |2005| ) and to evaluate errors in parameter estima 



tion by Monte-Carlo methods. Current drawbacks are the lack of 
proper halo-halo correlations and of large-scale velocity flows, 
and the absence of halo substructure for the resolved clusters. 
Note also that the approach based on a cluster mass function 
does not include the faint SZ emission of the filamentary struc- 
tures of the cosmic web. 

Comparison of cluster counts 



Figure 21 shows the cumulative counts for clusters of mass 
M500 > 10'** Mq, for a subset of the mass functions available 
to simulate the SZ effect and for different values of erg. Counts 
are quite different between the various options. The dependence 
of counts on erg is important and its precise value still uncertain, 
with a typical allowed range extending from about 0.81 as ob- 



tained from primary CMB constraints ( Komatsu et al. 



about 0.87 as obtained from weak lensing (Massey et a 



2 0TT|)to 
2007] ). 



Moreover, the dependence on the mass function is weaker but 
not negligible: for example, for the same set of WMAP 7-year 
cosmological parameters, the Jenkins and Sheth-Tormen mass 
functions predict about 50% more cluster above M500 > 10'"* Mq 
than the Tinker mass function. The sky model software can be 
used to produce maps with different mass functions and values 
of o-g. 

Constrained SZ simulations 

The simulation of cluster catalogues on the basis of a mass func- 
tion can be refined to constrain the catalogue to match available 
observational data. We then replace a subset of clusters of the 
randomly drawn catalogue by the clusters of the ROSAT all-sky 
survey and SDSS. The subset of replaced clusters in our cata- 
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logue is then chosen to match best the redshift-mass distribution 
of the MCXC and MaxBCG catalogues. The resulting maps can 
be computed quickly and are fully consistent with the X-ray and 
optical characteristics of the known clusters. 

5.1.3. SZ emission from large-scale structure simulations 

A complete and accurate simulation of the SZ effect requires 
a simulation of the distribution of hot electron gas in the uni- 
verse. As full scale simulations, including all the physics of the 
baryons in any given cosmological scenario, are not available at 
present, we adopt a simplified approach which relies both on A^- 
body+hydrodynamical simulations of the distribution of baryons 
for redshifts z < 0.025 (the local universe), and on pure A^-body 
simulations of dark matter structures in a Hubble volume. 

The local simulation of baryonic matter is obtained from the 
A^-body-i-hydrodynamical simulation of Dolag et al. (2005 ). The 
distribution is constrained by the observed IRAS 1 .2-Jy all-sky 
redshift survey, so that the positions and masses of the closest 
galaxy clusters are constrained to match quite closely those of 
real clusters. The distribution and peculiar motion of the bary- 
onic gas as obtained in the simulation permits us to simulate 
directly the corresponding thermal and kinetic SZ effect. 

The distribution of matter at higher redshift is obtained start- 
ing from a pure A^-body simulation of a Hubble-volume (a cube 
of 1 Gpc size). The SZ map for this latter simulation is obtained 
by putting, in the position of each dark matter concentration in 
the Hubble-volume simulation, a suitable cluster extracted from 
a template obtained using a different A'-body+hydrodynamical 
simulation of a smaller volume (Schafer et al.l 



of An was used for redshift radii up to z = 0.58. For redshifts 
exceeding z - 0.58, octant data sets, spanning a solid angle of 
7r/2, were replicated by rotation in order to cover the full sphere. 
Note that this generates fake replications of the farthest simu- 
lated clusters, the effect of which can be seen by careful inspec- 



Merging the two simulations, thermal and kinetic SZ maps 
are obtained by co-adding the contributions of the high redshift 
and low redshift structures. It is worth noting that while the 
higher redshift map contains only clusters, the hydrodynamical 
map used to describe the SZ emission from the local universe 
features the SZ from all the filamentary structure of the cosmic 
web. One pre-generated SZ simulation obtained in this way is 
currently available in the model. 

The A^-body-Hhydrodynamical SZ simulation has several ad- 
vantages. First of all, it is the only model presently implemented 
in the code that properly takes into account halo-halo correla- 
tions, models halo substructure and adiabatic gas physics, and 
provides cluster peculiar velocities in agreement with what ex- 
pected from the actual density contrast and its evolution with 
time. The simulation also models properly the filamentary struc- 
ture of the cosmic web for the low redshift universe. It also 
matches the observed low redshift matter distribution. The main 
drawbacks is the lack of flexibility: one single simulation is 
available, generated for fixed cosmological parameters {h - 0.7, 

= 0.3, Qa = 0.7, Ob = 0.04, trg = 0.9, «s = 1). In addition, 
the SZ map is generated on the basis of the extraction of a num- 
ber of halos from the output of the LSS simulation. Hence, the 
mass range of the clusters included in the map is limited by the 
cluster templates used: many low-mass systems, and even a few 
high-mass ones, may be missing in the final map. Finally, the 
cluster gas properties and scaling relations strongly depend on 
the physics included in the A^-body+hydrodynamical code used 
to generate the cluster templates. If some of these properties are 
not in complete agreement with X-ray cluster observations, this 
approach does not easily allow to modify them to correct for this 
effect. 

In order to cover redshifts out to the anticipated limit for 
Planck, several light-cone outputs were combined in the Hubble 
volume simulation: First, a sphere covering the full solid angle 



tion of the maps displayed in the left column of Fig. 22 



5.1 .4. SZ from LSS simulations and high z cluster counts 

In the previous approach, full hydrodynamical simulations at 
low redshift are complemented with a single high-redshift sim- 
ulation of large-scale structures. An alternate solution consists 
in merging a low redshift full hydrodynamical simulation, con- 
taining the constrained local SZ map, with a high redshift model 
based on cluster number counts. 

In this model, the low redshift SZ signal is simulated in con- 
centric layers around the observer using snapshots from full hy- 
drodynamical simulations up to z ^ 0.25. Seven layers, each 
with a comoving thickness of 100 /z"' Mpc, are used to generate 
the required volume. The thermal and kinetic SZ signals are then 
integrated using the method in da Silva et al. (2000 2001 1 and 
the HEALPix sky pixelisation scheme. For the innermost layer 
we used the local constrained simulation in Dolag et al. ( 2005 [ l, 
i.e., the baryon distribution and the resulting local constrained 
SZ map used in Section 5.1.3 For the remaining layers up to 
z - 0.25 we use the ACDM simulation inl De Boni et al.| ( |201()l l. 
Both hydrodynamical simulations feature full baryonic gas dy- 
namics and state-of-the-art modelling for gas cooling, heating by 
UV, star formation and feedback processes. 

We then draw at random additional clusters modelled as in 



Section 5.1.2 for z > 0.25. As in the previous approach, the sim- 
ulation is generated for fixed cosmological parameters (h = 0.70, 
= 0.27, Qa = 0.73, Ob = 0.044, erg = 0.78, «s = 0.95). The 
main advantage of this approach with respect to the previous 
one is that the simulation features a model for the local SZ dif- 
fuse component and filamentary structures up to z = 0.25. Note 
that the cosmological parameters used are different from those 
of Section |5T3] 

5.1.5. Polarised SZ effect 

The polarised SZ effect arises from local quadrupoles in the 
incoming radiation pattern as seen from the cluster (Audit & 
Simmons 1999; Sazonov & Sunyaev[ [T999| |Liu, da SilvaTfe 
Aghanim 2005). Such quadrupoles can be either cosmological 
(the CMB quadrupole as seen from the location of the cluster), or 
kinetic (due to a relativistic effect of the transverse motion of the 
cluster). This effect is neglected in the current implementation. 



5.1 .6. Limitations of the SZ models 

The SZ maps derived from A^-body+hydrodynamical simula- 
tions are stored in HEALPix format, with A^side = 1024 (pixel 
size of 3.44') and A^side = 2048 (pixel size of 1.72') respec- 
tively. Running the PSM at smaller pixel size would not upgrade 
the resolution of the maps, i.e., in spite of the larger number of 
pixels, no cluster substructures at scales below 1.72' would be 
present. 

The semi-analytic model does not suffer from this limitation, 
as clusters are modelled using analytic profiles which are scaled 
and normalized using scaling laws between flux, size and red- 
shift, mass. This model, however, does not include cluster sub- 
structure, nor scatter in the scaling relations and profiles. Hence, 
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Fig. 22. Maps of thermal and kinetic SZ effect from the Hubble volume simulation (left column), from the local universe (middle 
column), and total thermal and kinetic SZ effects from both simulations together, smoothed here to a sky resolution of 5 arc-minutes 
(right column). All the displayed maps are small 5° x 5° patches centred on the North Galactic Pole. The original full sky maps are 
at HEALPix A^side - 2048 for the Hubble volume, and at A^side = 1024 for the local hydrodynamical simulations. The contribution 
from the local universe comprises a large cluster at sky coordinates close to those of the Coma cluster, which dominates the total 
SZ emission in the maps displayed here. 



two distinct clusters that have the same mass and redshift would 
be strictly identical on the modelled SZ map. In addition, the dif- 
fuse SZ effect is not modelled in the current version, as the SZ 
effect is generated by spherically symmetric clusters distributed 
on the sky with uniform probability and no correlation. 



Figure 23 illustrates the large fluctuations these different op- 



tions yield in terms of the thermal and kinetic SZ power spectra. 
The semi-analytic model can be adjusted to fit current measure- 
ments at high (, adjusting the value of cr% and the Y-M scaling 
law, but has significantly less power on large scales than simu- 
lated SZ maps based on hydrodynamical simulations at low red- 
shift. 



The current implementation of the thermal SZ effect is based 
on maps of the Compton y parameter. The relativistic SZ effect 
is included only at some level for the models that are based on 
the generation of a cluster catalogue. 



5.2. Radio sources 

Radio sources in our model are mainly modelled on the basis of 
surveys of radio sources at frequencies ranging from 0.85 GHz 
to 4.85 GHz. A summary of these surveys is given in Table |2] 
Figure |24]illustrates their sky coverage. For the present purpose 
it is useful to distinguish four cases: 

1. green points are sources with fluxes measured at both ap- 
proximately 1 GHz (NVSS or SUMSS) and 4.85 GHz (GB6 
orPMN): 109,152 sources; 

2. blue (yellow) points are sources present only in the NVSS 
(SUMSS) and thus with fluxes only at 1.4 (0.843) GHz: 
1,813,551 (158,284) sources; 

3. red points are sources with only PMN (4.85 GHz) fluxes: 
8656 sources; 

4. the small white regions are those not covered by any survey. 

Extrapolation of the flux of each source to all useful frequencies 
requires the knowledge of the spectral behaviour of the sources 
which, for any individual object frequently is quite complex (see 
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Table 2. Summary of the large-area surveys of point sources used in this work. 



Frequency 


Catalogue 


J lim(rnjy) 


Dec ran: 




Angular resolution 






GB6 


18 





+75 


3.5' 


Gregory et al. (1996) 




PMNE 


40 


-9.5 - 


+ 10 


4.2' 


Griffith et al. (1995) 


4.85 GHz 


PMNT 


42 


-29 


-9.5 


4.2' 


Griffith et al. (1994) 




PMNZ 


72 


-37 


-29 


4.2' 


Wright etal. (1996) 




PMNS 


20 


-87.5 - 


-37 


4.2' 


Wright etal. (1994) 


1.4 GHz 


NVSS 


2.5 


-40 


+90 


45" 


Condon etal. (1998) 


0.843 GHz 


SUMSS 


18 


-50 


-30 


45" cosec(|(5|) 


Mauch et al. (2003) 
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-50 


45" cosec(|(5|) 
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Fig. 23. Power spectra De - (((+ l)C[/2n of thermal and ki- 
netic SZ effect, as modelled by the PSM for various options. 
Black (resp. green) curves give the spectrum of the tSZ (resp. 
kSZ) maps as modelled on the basis of a population of clusters 
as described in 5.1.2 Upper curves are obtained using the mass 
function of Tinker et al.| ( |2008" i with WMAP 7-year cosmologi- 



cal parameters and the Y-M scaling law of Arnaud et al. ( 2010[ l. 
Lower curves are obtained when the value of erg is set to 0.75 
instead of 0.809, and assuming an hydrostatic bias of 15% in 
the Y-M scaling law, i.e., the Y parameter is 15% lower than 
in Arnaud et al. ( 2010 1 for a given cluster mass. Red (resp. blue) 
curves are tSZ (resp. kSZ) power spectra for the model described 
in 5.1.4 where the high redshift part is generated with the same 
two modeling alternatives (upper curves for cr^ = 0.809 and no 
hydrostatic bias, lower curves for erg = 0.75 and 15% hydro- 
static bias). As the low-redshift part of the SZ maps uses maps 
at HEALPix A^side - 1024, power spectra for this case are com- 
puted only up to ^ = 2000. The red diamond at { - 3000 is 
the tSZ measur ement of D'-^^ - 3.65 + 0.69 /l/K^ obtained re- 
cently by SPT ( [Reichardt et al.| 
measurements from ~ 



Lueker et a 



201 l|l, improving on previous 
' pOTO| ), |Dunkley et al.| ( |201 1| , 



and Shirokoff et al. ( 2011 1. The blue arrow is the corresponding 
upper Umit for kSZ. 



e.g. Sadler et al. 2006). For the present work, we resort to a 
power law approximation of the spectra of the sources, of the 
form 5v oc v^". Spectral indices, however, are different in a few 
pre-selected intervals of electromagnetic frequency. 

For sources of the group (1) above, which have measure- 
ments at two frequencies, the individual spectral indices a below 
5 GHz can be directly estimated as 




150 2C0 



a = -log(54.85GHz/5low)/log(4.85/viow), 



(51) 



Fig. 24. Sky coverage of the surveys listed in Tabled in Galactic 
coordinates. Green points: sources present in both ^ 1 GHz 
(NVSS or SUMSS) and 4.85 GHz (GB6 or PMN) catalogues; 
blue points: sources in the NVSS catalogue only; yellow points: 
sources in the SUMSS catalogue only; red points: sources in the 
PMN catalogue only; white regions: not covered by any survey. 



where vio„ is either 1 .4 GHz, if we are in the region covered by 
the NVSS, or 0.843 GHz if we are in the lower declination re- 
gion covered by the SUMSS. However, the calculation is not as 
straightforward as it may appear, because the surveys at differ- 
ent frequencies have different resolutions, implying that a single 
source in a low resolution catalogue can correspond to multi- 
ple sources at higher resolution. The spectral index estimates 
were carried out by degrading the higher resolution survey to 
the resolution of the other In practice, whenever the higher res- 
olution (NVSS or SUMSS) catalogue contains more than one 
source within a resolution element of the 4.85 GHz one, we have 
summed the NVSS or SUMSS fluxes, weighted with a Gaussian 
beam centred on the nominal position of the 4.85 GHz source, 
and with FWHM equal to the resolution of the 4.85 GHz survey. 

On the other hand, the low frequency surveys, and espe- 
cially the NVSS, are much deeper than the 4.85 GHz ones, which 
also have inhomogeneous depths. Simply summing all the lower 
frequency sources within a resolution element has the effect 
of including a variable fraction of the background, unresolved 
at 4.85 GHz, thus biasing the spectral index estimates. To cor- 
rect for this, we have selected 159,195 control fields, free of 
4.85 GHz sources, and computed the average flux of NVSS or 
SUMSS sources within a 4.85 GHz beam pointing on the field 
centre, again taking into account the 4.85 GHz response func- 
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Fig. 25. Modelled source number counts at 5 and 20 GHz for 
one sky realisation, normalised to AA^o - S(]yT^'^, compared 
with models and observational data. Data at 5 GHz are from 
Kellermann et al.| (1986), F omalont et al.| ( |199l"] l, and |Haarsina| 
et al. (2000 1. Data in the 20 GHz panel are from the 9C survey 
( |Waldram etaL|[2003) at 15 GHz and from the ATCA survey at 
8 GHz ( |Ricci et al.||2004a| l; no correction for the difference in 
frequency was applied. 



Fig. 26. Comparison of modelled radio source counts (red 
points) and the Planck radio counts (blue points, IPlanck 



[Collaboration Xni]|201 l| l. Also shown are: the counts estimated 
at 31 GHz from DASI ( grey dashed box [Kovac etaL||2002l ) and 
PACO (grey diamonds, Bonavera et al. 201 l| l, at~33GHz from 



the VSA data (gre y box,|C leary et al.||2b05| l7an d the SPT (grey 



squares, [Vieira et al . 2010 ) and ACT (grey stars. Marriage et al. 
|201 lb| l counts of radio sources. 



tion. The average fluxes of control fields, ranging from 1.67 to 
3.16 mJy for the different combinations of low and high fre- 
quency catalogues (NVSS/GB6, NVSS/PMN, SUMSS/PMN), 
have been subtracted from the summed NVSS or SUMSS fluxes 
associated with 4.85 GHz sources. In this way we obtained spec- 
tral indices from 1 to 5 GHz for a combination of complete 
5 GHz selected samples with somewhat different depths, sum- 
ming up to 109,152 sources over about 95% of the sky. 

Although this is the best we can do with the available data, it 
is clear that the derived individual spectral indices are uncertain, 
due to a combination of several factors: measurement errors, un- 
certainties associated to the corrections applied, that are of sta- 
tistical nature, and variability (the surveys have been carried out 
at different epochs). As a result, the absolute values of several 
individual spectral index estimates turned out to be unrealisti- 
cally large, and the global distribution was found to be much 
broader than indicated by accurate studies of smaller samples. 
Furthermore, we obtain an estimate of the 5 GHz counts exceed- 
ing by an average factor of roughly 1.8 those directly observed, 
again indicating that the spectral index distribution is spuriously 
broadened. Therefore, we use these spectral index estimates only 
to assign the sources either to the steep- or to the flat-spectrum 
class (the boundary value being a - 0.5) and to determine the 
mean spectral index of each population. We find (asteep) - 1-18, 
(ttflat) = 0.16. The spectral index distributions are approximated 
by Gaussians with variances crjteep.flat = 0.3, consistent with the 
results by Ricci et al.| (p006|. We then extrapolate the 5 GHz 
fluxes to 20 GHz by assigning to each source a spectral index 
randomly drawn from the Gaussian distribution for its class. 

Sources with flux measurements at a single frequency have 
been randomly assigned to either the steep or to the flat-spectrum 
class in the proportions observationally determined by Fo malont^ 
et al. ( 1991| l for various flux intervals, and assigned a spectral 
index randomly drawn from the corresponding distribution. The 



holes in the sky coverage have been filled by randomly copy- 
ing sources from other regions in proportion to the surface den- 
sity appropriate for each flux interval. The same procedure was 
adopted to add fainter sources in the regions where the existing 
surveys are shallower, until a coverage down to at least ^ 20 mJy 
at 5 GHz over the full sky was achieved. We have checked that 
still fainter sources do not appreciably contribute to fluctuations 
in Planck channels for detection limits in the estimated range 
(approximately 200 to 500 mJy; Lopez-Caniego et al. 2006), as 
expected since fluctuations are dominated by sources just be- 
low the detection limit. This check is caiTied out computing the 
power spectra of fluctuations due to sources below such limits 
in regions covered by the NVSS (the deepest survey) and in re- 
gions covered to shallower limits: the results are indistinguish- 
able. This, however, holds for simulations of Planck observa- 
tions. For other experimental setups, this check should be re- 
peated considering a correspondingly different combination of 
resolution and sensitivity. 

The regions less covered by the surveys used in the model, 
where the fraction of simulated sources is larger, are mostly 
around the Galactic plane, where they have a small effect com- 
pared to free-free and synchrotron emissions. At Galactic lati- 
tudes \b\ > 10° the fraction of real sources (at least as far as po- 
sitions are concerned) is about 97%; over the full sky it is about 
95%. Therefore we expect that the simulated maps faithfully re- 
flect also the clustering properties of radio sources. 

In Figure |25] the source counts at 5 and 20 GHz obtained 
from our map are compared with observed counts, with the 



model by |Toffolatti et al. ( 1998 1, and with an updated version of 
the model by lde Zotti et al.| ( |2005 ), allowing for a high-redshift 
decline of the space density of both flat-spectrum quasars (FSQs) 
and steep-spectrum sources (not only for FSQs as in the origi- 
nal model). The model adopts luminosity functions (in units of 
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Fig. 27. Emission law for a selected sample of typical radio 
sources produced in a realisation of the sky emission with our 
model. 
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(52) 



and lets those of steep-spectrum sources and of FSQs evolve in 
luminosity as 



i*,FSRQ(z) = L.(0)10' 



,*:cvZ(2z,„p-z) 



(53) 



while for BL Lac objects a simpler evolutionary law is used 
L,(z) = L,(0)exp[^evT(z)], 



where t(z) is the look-back time in units of the Hubble time, //q 
and where fegv and Ztop parametrize the luminosity evolution. The 
new values of the parameters are given in Table[3] Figure 26 sim- 
ilarly shows a comparison of the radio source number count in 
the present model with predicted number counts from the Degree 
Angular Scale Interferometer (DASI) and the Very Small Array 
(VS A), and observed number counts from Planck, Planck ACTA 
Coeval Observations (PACO), the South Pole Telescope (SPT) 
and the Atacama Cosmology Telescope (ACT). 

Note that we do not distinguish, for the moment. Galactic ra- 
dio sources from extragalactic ones. While comparisons of num- 
ber counts are made with cosmological evolution models for ex- 
tragalactic radio sources. Galactic sources do affect only number 
counts at low galactic latitudes, typically \b\ < 5°. 

Due to the complex spectral shape of radio sources, the 
power-law approximation holds only for a limited frequency 
range. To extrapolate the fluxes beyond 20 GHz we use the mul- 
tifrequency WMAP data (Bennett et al. 2003) to derive the dis- 
tributions of differences, Sa, between spectral indices above and 
below that frequency. Such distributions can be approximated 
by Gaussians with mean 0.35 and dispersion 0.3. To each source 
we associate a spectral index change drawn at random from the 
distribution. Finally, for a small number of sources with exceed- 
ingly low negative values of a, yet another spectral index is as- 
signed at frequencies higher than 100 GHz. In summary, each 
radio source in the model is simulated on the basis of four dif- 
ferent power laws: One at frequencies below 5 GHz, another one 
between 5 and 20 GHz, another one between 20 and 100, and a 
last one at v > 100 GHz (which is, for most sources, the same as 
between 20 and 100 GHz). 

We also produce maps and catalogues of polarised radio 
source emission attributing to each source a polarisation degree 



Fig. 28. Similar to 26 but with additional green points showing 
the counts of WMAP radio sources as represented in our model. 
The solid black curve shows the total number counts of extra- 
galactic radio sources as predicted by the updated model of de 
Zotti et al. (2005). 



randomly drawn from the observed distributions for flat- and 

2004b| l, and a 



steep-spectrum sources at 20 GHz (Ricci et al 



(54) polarisation angle randomly drawn from a uniform distribution. 



5.2.1 . The special case of WMAP sources 



Sources available in the WMAP Point Source catalogue (Gold et 
|aL]|201 1| ) have been included in the model as follows: 

Given a WMAP source, we search for a match in the avail- 
able radio sources in the catalogue built from low frequency sur- 
veys. The identification is made with the brightest radio source 
at 4.85 GHz within 11' from the WMAP source location. This 
method provides radio counterparts for all the WMAP sources. 
The value of 11' coincides with the one used by WMAP authors 
for source identification at 5 GHz. 

The matching radio sources are removed from the radio 
source catalogue described in Section 



5.2 and are included in 



a separate catalogue which merges measured low frequency and 
WMAP fluxes. 

These WMAP sources are simulated with seven power laws: 
The first at frequencies below 5 GHz, and then between 5 and 
23 GHz {WMAP K band), between 23 and 33 GHz (Ka band), 
between 33 and 41 GHz (Q band), 41 and 61 GHz (V band), 61 
and 94 GHz (W band), and the final emission law for frequencies 
over 94 GHz. In this last frequency range, the spectral index is 
taken to be the same as between 61 and 94 GHz, unless it is 
positive (flux increasing with frequency) in which case it is set 
to 0. 

Note that most of WMAP radio sources are highly variable 
( Gonzalez-Nuevo et al. 2008| l. Including them in the simulations 
hence makes sense for simulating WMAP sky observations, but 
not for all types of simulations. 

If the simulation is polarised, then the WMAP sources polar- 
isation is set at the (simulated) value of their radio counterpart. 
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Table 3. Best-fit values of the parameters of the evolutionary models for canonical radio sources. 



Source type 






luminosity function 




evolution 




log;;o(Mpc 


a 


b log L, (ergs" 


'Hz ' at5GHz,z = 0) 


^ev ^top 


FSQ 


-8.989 


0.658 


2.938 


34.043 


0.224 2.254 


BL Lac 


-7.956 


0.975 


1.264 


32.831 


1.341 


Steep 


-7.389 


0.729 


2.770 


33.177 


0.262 2.390 



5.2.2. Limitations of tine radio source model 

Extrapolations of radio source counts to frequencies > 30 GHz 
show evidence of substantial incompleteness below about 0. 1 Jy 
(see Fig. |28i. This is not a serious problem for WMAP or 
Planck-related simulations because the fluctuations due to un- 
resolved radio sources are dominated by sources brighter than 
this limit. On the other hand, the current version of the PSM un- 
derestimates the amplitude of radio shot noise for experiments 
with higher spatial resolution and correspondingly lower con- 
fusion noise levels. Note, however, that at frequencies greater 
than 200 GHz, the radio shot noise becomes quickly negligible 
compared to that due to far-IR galaxies. In the next release of 
the PSM this problem will be cured by taking into account the 
data from the AT20G survey ( [Hancock et al. 2011; Murphy et 
|aLjp010 ), that allows a much better determination of the flux- 
density dependent distribution of spectral indices up to 20 GHz, 
as well as the high frequency counts and spectral information 
from the SPT survey at 150 and 220 GHz ([Vieira et al | [20T0l l 
and from the ACT survey at 148 GHz ( [Marriage et al.||2011b| l. 



Further information on high-frequency source spectra is pro- 
vided by the quasi-simultaneous multifrequency observations re- 
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ported by Bonavera et al. (201 1 1; Massardi et al. (2011 1; Planck 



; [Planck CoUaboration XV[ ( [201 l| l; 



5.3. Far-Infrared Sources 
5.3.1. IRAS sources 

The software uses a compilation including all FIR sources taken 



from the IRAS Point Source Catalog (PSC, Beichman et al. 



1988|) and the Faint Source Catalog (FSC, [Moshir, Kopman, & 
Conrow 1992| l. Sources identified as ultra-compact Hii regions 
are removed from the catalogue, and given a special treatment 
(see Section [54] l. 

For the remaining sources, fluxes are extrapolated to Planck 
frequencies, adopting modified blackbody spectra i/B(v, T), 
B(v, T) being the blackbody function. For sources detected at 
only one IRAS frequency, the values of b and T are taken to be 
those of the average spectral energy distribution of the sample 
by [Dunne et"ni| ( [2g00| ), i.e., b = 1.3 and T = 35K. If the source 
is detected at 60 and 100 fim, then b = 1.3 is stiU assumed but 
the temperature is obtained from a fit to the data. 

Since the PSC does not contain objects where the confusion 
from Galactic sources is high, and thus does not penetrate the 
Galactic centre well, and since the FSC is restricted to regions 
away from the Galactic plane, the source density is a function 
of Galactic latitude. For each PSM simulation, we add randomly 
distributed sources until the mean surface density as a function 
of flux matched everywhere the mean of well covered regions 
down to S 857GHz ~ 80 mJy. The coverage gaps of the catalogue 
(the IRAS survey missed about 4% of the sky) are filled by the 
same procedure. 



To each source we assign a polarisation fraction drawn from 
a distribution with one degree of freedom, the global level of 
which can be adjusted as an input parameter (and is 1% by de- 
fault). The polarisation angle is randomly drawn from a uniform 
distribution. 

In the near future, the catalogue of observed far-infrared 
sources implemented in the model will be complemented with 
the Planck Early Release Compact Source Catalogue (ERCSC, 
[Planck Collaboration Vn[[20n] l. 



5.3.2. Cosmic Infrared Background anisotropies 

An important, and possibly dominant, contribution to sub- 
millimetre small-scale anisotropies comes from galaxies se- 
lected by SCUBA and MAMBO sui-veys (see, e.g., [Scott, 



Dunlop, & Serjeant 2006 ; Coppin et al. 2006 ), that are probably 
strongly clustered (Negrello et al. , 2004, 2007). These galaxies 
are interpreted as massive proto-spheroidal galaxies in the pro- 
cess of forming most of their stars in a gigantic starburst. We 
adopt the counts predicted by the model of |Lapi et al. ( 2006| l, 
which successfully accounts for a broad variety of data includ- 
ing the SCUBA and MAMBO counts and the preliminary red- 



shift distribution ( jChapman et al. 2005 1. 

The clustering properties of these sources are modelled as 
in [Negrello et al. (2004), using their more physical model 2. 
The simulation of their spatial distribution is produced using 
the method by Gonzalez-Nuevo, Toffolatti, & Argiieso ( 2005| l. 
Once a map of the source distribution is obtained at the ref- 
erence frequency Vref, extrapolations to any other frequency v, 
are obtained via the flux-dependent effective spectral indices 
a - -log(5ref/5,)/log(vi-ef/v,), where S, is defined by n(> 
5,;v,) - n(> 5,ef;vref). The spectral indices are computed in 
logarithmic steps of Alog(5ref) = 0.1. We check that the counts 
computed from the extrapolated maps accurately match those 
yielded, at each frequency, by the model. 

As first pointed out by Blain[ ( [T996 |, the combination of the 
extreme steepness of the counts determined by SCUBA surveys 
and of the relatively large tensing optical depth corresponding 
to the substantial redshifts of these sources maximises the frac- 
tion of strongly tensed sources at (sub)-mm wavelengths. We 
include such sources in our simulation by randomly distribut- 
ing them with flux-dependent areat densities given by [Perrotta et[ 
al. (2003). The frequency extrapolations are made via the spec- 
tral indices obtained in the same way as for the proto-spheroidal 
galaxies. 



Figure 29 shows the power spectrum of cosmic infrared 
background anisotropies as computed on a simulated sky, and 
as measured by Planck (Planck Collaboration XVIII |201 l| l, and 
ACT (Dunkley et al. 201 l| l. The^ agreement between the model 
and the measured power spectra is reasonable, but not perfect, 
and discrepancies will be addressed with improved modeling in 
future releases, in the light of recent advances on number counts 
and on correlations from observations with the BLAST balloon- 
borne experiment ( [Patanchon et al.[ [2009[ [Viero et"aL| [2009| l, 
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Fig. 29. Cosmic infrared background power spectrum. Solid 
lines are obtained from a simulation. Data points are from Planck 
observations ( jPlanck Collaboration XVIII] |201 l| l. Dashed lines 
at high multipoles for 143 and 217 GHz are from the Dunkley 
[etaLlpOTT] ) best-fit model for IR sources (at 148 and 218 GHz, 
extracted from Fig. 2 of their paper). The CMB power spectrum 
(in black) is a theoretical model fitting WMAP observations. 



SPT (IHall et"aLl[20T0| , Herschel (Maddox et al.| [20T0l[Cooray 



etal. 



2010 1, and from a combined analysis of ACT and BLAST 



( Hajian et al. 2012 1. Also, the correlation between the CIB maps 
across the frequency range is almost unity in the current model, 
a feature that will have to be improved in the future using a more 
detailed model constrained by upcoming additional Planck ob- 
servations. 

Cosmic infrared background anisotropics as implemented in 
the current model are assumed unpolarised. 

5.3.3. Limitations of fine far-infrared source model 

The main uncertainties in the model for relatively nearby far- 
infrared galaxies arise from the extrapolations of IRAS flux 
densities. The Planck data ( [Planck Collaboration XVI| |2011| l 
show evidence for colder dust than has previously been found 
in external galaxies so that our model likely underestimates the 
(sub)-mm flux densities of at least a fraction of IRAS galaxies 
and, consequently, the source counts. Moreover, the CO emis- 
sion lines are not included in the simulations. The CO(J=l— >0), 
CO(J=2^1), CO(J=3^2), fliat, for low-z galaxies, fafl within 
the Planck 100, 217 and 353 GHz passbands, respectively, may 
contribute significantly to the observed fluxes, while higher or- 
der CO transitions are probably negligible. A careful investiga- 
tion of the possible contamination by these lines is in progress 
in preparation for future releases of the PSM. 

A wealth of data on faint far-IR sources, those that dominate 
the small-scale fluctuation at frequencies higher than 200 GHz, 
have become available after the extragalactic source model was 



worked out. These include source cou nts ^Clements et al. 



[Oliver et al.| [20T0l [Blhermin et al.| [2010t [Glenn eFal 



2010 



2010 



and estimates of the power spectrum of the cosmic infrared back- 
ground anisotropics (Planck Collaboration VIII 201 1[ Amblard 



et al. 201 1 1. The models implemented in the PSM fare reason 



of dusty galaxies, coupled with the fast cosmological evolution, 
make both the faint counts and the CIB power spectrum exceed- 
ingly sensitive to the detailed spectral energy distribution (SED) 
of model galaxies. For example, the amplitude of the (sub-)mm 
power spectrum scales as the frequency to a power of 7-8 so that 
relatively small differences in the SED are enough to yield large 
discrepancies between the model and the data. Much improved 
models, successfully reproducing the recent data, are now avail- 
able both for the epoch-dependent luminosity function (hence 
for the source counts; e.g. |Bethermin et al. 2011 Lapi et al. 



2011 1 and for the clustering properties (e.g. Xia et a 



2011 



Renin et al. 2012). These models will be exploited in next re- 



leases of flie PSM. 

5.4. Galactic point sources 

Most Galactic compact sources are currently treated in the model 
as part of the diffuse emission. IR sources present in the IRAS 
catalogue are treated as extragalactic sources (dusty galaxies). 
A special treatment has been given to simulate emission from 
ultra-compact H n (UCHii) regions, however 

A list of UCH II regions and UCH n region-candidates has 
been extracted from the IRAS Point Source Catalogue with 
fluxes from radio counterparts. A total of 864 sources were se- 
lected from IRAS PSC according to the following criteria. 



1 . IRAS colours as specified in [Kurtz, Churchwell, & Wood 
(jl994j). 

2. Fluxes in the IRAS channels 100 /urn, 60 fim and 25 i-im are 
not upper limits (i.e., they are at least of "medium" quality. 



according to definitions given in [Beichman et al. ( 1988) ). 

3. The flux at 100 yum was required to be greater than or equal 
to 100 Jy. 

4. The Galactic latitude is below 10°. 

This list includes all confirmed UCH n regions from the orig- 
inal sample of Kurtz, Churchwell, & Wood (1994) that have an 
IRAS counterpart (48 objects) and 38 UCHii regions out of 53 
of the sources in |Wood & Churchwell| ( |1989| l. 

Note that as some of these sources are present in the dust 
map used for modelling thermal dust emission, that map has 
been processed to subtract the UCH n regions, to avoid double- 
counting. 

Radio counterparts have been searched for in the GB6 and 
NVSS catalogues. The search radius is 1'; when more than one 
counterpart was found within the search radius (within that cat- 
alogue), the brightest was selected. None of these radio cata- 
logues covers the Southern hemisphere: search for radio coun- 
terparts in catalogues of the southern hemisphere is in progress. 

For all 864 sources, the fluxes in the useful frequency range 
are estimated by performing a modified blackbody fit to the 
IRAS fluxes at 100 and 60 fim. The emission is parametrised 
with the form 



S ^vf^ Biv) = v' 



2hv^ 



1 



exp{hv/kT) - 1 



(55) 



with jg ^ 1.5 (see e.g. IGear, Robson, & Griflin"] [1988) [Hoare, 



ably well with the new data, but the steep (sub-)mm spectrum 



[Roche, & Glencros s 1991 ; Maxl a et al.|[2001[ )r 

When a radio counterpart to the IRAS source is found, fluxes 
at low frequencies are corrected by adding to the modified black- 
body fit an estimated low-frequency flux extrapolating the flux 
at the frequency of the radio catalogue (where the counterpart 
was found) with a free-free spectral index. When a radio coun- 
terpart to the IRAS source is found in more than one of the three 
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radio catalogues priority was given to GB6, then to |Giveon et al.| 
( |2005l l and then NVSS. 

In case of an erroneous radio counterpart this procedure 
may lead to slightly overestimated fluxes at the lowest HFI fre- 
quencies, but this has little impact for frequencies greater then 
200 GHz, where we expect these sources to be very relevant for 
Planck. 

UCH II regions are free-free emitters and lack magnetic fields 
to cause grains to align. Hence, they are not intrinsically po- 
larised, and are thus modelled as unpolarised in the current ver- 
sion of the model. 

6. Conclusions 

We have developed a complete, flexible model of multicompo- 
nent sky emission, which can be used to predict or simulate 
astrophysical and cosmological signals at frequencies ranging 
from about 3 GHz to 3 THz. The model, which actually imple- 
ments several options for each of the components, has been de- 
veloped for testing component separation methods in prepara- 
tion for the analysis of Planck data, but has been used also for 
simulating data in the context of analyses of WMAP observa- 
tions, or for planning future missions such as COrE ( The COrE 
Collaboration et al. |201 l| l. Given the high sensitivity of these 
instruments, the quality of the reconstructed CMB temperature 
map depends strongly on our ability to remove the contamination 
by astrophysical foregrounds. It is therefore necessary to build 
simulations as close as possible to the complexity of the real sky 
over a large frequency range. We present here what has been 
achieved through a series of improvements accomplished over 
several years and the contributions of experts in different fields: 
CMB, Galactic emission and compact sources, extragalactic ra- 
dio and far-IR sources, Sunyaev-Zeldovich signals from clusters 
of galaxies and the intergalactic medium. Note however that the 
current version of the model cannot be expected to match per- 
fectly data sets that have not been used in the model, in particular 
those of the Planck mission and of other upcoming observations. 
Updates will be proposed as such data sets become available. 

The maps of each astrophysical component have been re- 
peatedly updated by requiring a close match with the constantly 
increasing amount of observational data. Fig. 30 shows the spec- 



trum of root mean square (rm.s.) fluctuations in simulated maps 
at WMAP frequencies for \b\ > 20° and 10° resolution. The var- 
ious diffuse components have different spectral characteristics 
while the total signal is a good match to the WMAP 7 -year data; 
the rm.s. values agree to within better than 5 percent. The sky 
emission maps observed with WMAP, smoothed to 1° resolu- 



tion, are compared in Fig. 31 to PSM predicted emission maps 
at the same frequencies and resolution. The agreement is excel- 
lent over most of the sky, except in the galactic ridge and around 
a few compact regions of strong emission, for which the model 
does not exactly reproduce the observations. Discrepancies are 
due to a mixture of uncertainties in the exact resolution of the 
model (and, in particular, a lack of resolution of the map of syn- 
chrotron spectral index), errors in extrapolation of the dust emis- 
sion, residuals of galactic emission in the predicted CMB map, 
and insufficient number of parameters used in the current PSM in 
regions of strong emission (in particular when emission comes 
from several distinct regions along the line of sight). The model 
will be refined in future version, in particular using Planck ob- 
servations. 

The CMB temperature and polarisation maps, currently 
based on WMAP observations (maps and best-fit cosmological 
model), are complemented with simulations to allow for the 
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Fig. 30. Temperature r.m.s. fluctuations at WMAP frequencies 
for \b\ > 20° at a resolution of 10°. The symbols represent the 
fluctuations in the various diffuse components of the sky model, 
the total simulated fluctuations, and WMAP 7-year maps. The 
total signal is a good match to the WMAP data. 



weak gravitational lensing by gradients in the large-scale gravi- 
tational field. Weak gravitational lensing has been dealt with in 
the Born approximation. Simulations of CMB temperature maps 
with primordial non-Gaussianity of local type are also accessible 
via our sky model. 

The Galactic diffuse emission model includes five compo- 
nents: synchrotron, free-free, spinning dust and thermal dust ra- 
diation, and '^CO (J=l-0), (J=2-l), and (J=3-2) molecular fines 
at 115.27, 230.54, and 345.80 GHz, respectively. For each com- 
ponent, alternative models are available but we propose a com- 
bination of models that reproduces best the available data. In this 
default model the synchrotron emission is based on the 408 MHz 
all-sky map by Haslam et al. (T982| extrapolated in frequency 
exploiting the spectral mdex map by Miville-Deschenes et al. 
(2008; model 4). The free-free template is obtained from the 
WMAP MEM map, complemented with the Ho- all-sky template 
by Finkbeiner (2003) corrected for extinction as in Dickinson et 
al. (2003). For thermal dust we adopted model 7 of Finkbeiner 
et al. (1999), which features spectral variations across the sky. 
For spinning dust we adopted the template produced by Miville- 
Deschenes et al. (2008) from an analysis of WMAP data. The CO 
emission map is based on the '^CO(l-O) survey by Dame et al. 
(2001). Standard intensity ratios with the (J=l-0) line have been 
used for the (J=2-l) and (J=3-2) transitions. Since all these tem- 
plates have a resolution insufficient for our purposes, the code 
allows the possibility to add small-scale fluctuations following 
Miville-Deschenes et al. (2007). 

In our model as currently implemented, the only diffuse 
emissions that are polarised (apart from the CMB) are syn- 
chrotron and thermal dust. For synchrotron we rely on the 
WMAP 23 GHz polarisation maps, extrapolated in frequency us- 
ing the same spectral index map used for intensity. Modeling 
dust polarisation is made difficult because of the paucity of 
data. In our model the dust Q and U maps are built assuming 
a constant intrinsic polarisation fraction, geometrical depolari- 
sation and polarisation angle maps constructed using a Galactic 
magnetic field model for scales larger than 20° and constraints 
from 23 GHz WMAP polarisation data to reproduce the projected 
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Fig. 31. Comparison of sky emission as observed by WMAP (7-year data) and as predicted by the PSM in the same frequency bands, 
at a resolution of 1°. For each frequency channel, the color scale is the same for WMAP (left column) and PSM prediction (middle 
column). An histogram equalised color scale is used for the K, Ka, and Q channels, and a linear scale for the V and W channels. 
Maps are saturated to highlight common features away from the galactic plane. Maps of difference between PSM prediction and 
WMAP observation are displayed in the right column (note that the color scale is different from that used to display the K, Ka and 
Q maps), highlighting discrepancies in the galactic plane specifically and at the location of a few regions of compact emission. The 
agreement is excellent over most of the sky away from the galactic ridge and a few compact regions. 



structure at smaller scales due to the turbulent magnetic field of 
the ISM. 

While generally compact Galactic sources are part of the 
diffuse emission map and Galactic point sources are treated in 
the same way as extragalactic point sources (catalogues do not 
distinguish between Galactic and extragalactic point sources), 
ultra-compact H ii regions have been included in the model as 



a separate population. They have been selected from the IRAS 
catalogue and the extrapolation in frequency of their flux densi- 
ties is made adopting a grey-body spectrum. Radio counterparts 
have been searched for in the NVSS and GB6 catalogues. 

The Sunyaev-Zeldovich effect from galaxy clusters is sim- 
ulated in two ways. One can start from the epoch-dependent 
cluster mass function, for the preferred choice of cosmological 
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parameters, plus some recipes to model the density and tem- 
perature distribution of hot electrons. A fraction of simulated 
clusters can be replaced by real clusters drawn from the ROSAT 
all-sky or the SDSS catalogue. Alternatively, we can resort to 
A^-body-hhydrodynamical simulations that also contain the SZ 
emission from the web of intergalactic hot gas. 

Most radio sources in the model are real, taken from the rela- 
tively deep all-sky low-frequency catalogues. Each source has its 
own spectral properties, either determined directly from multi- 
frequency data or randomly extracted from the observed distri- 
butions of spectral indices. Extrapolations to higher frequencies 
are made using different spectral indices for different frequency 
intervals, in order to ensure consistency with multi-frequency 
source counts. This approach turned out to be remarkably suc- 
cessful in predicting the counts that have been later measured 
by Planck. The polarisation degree of each source is randomly 
drawn from the distributions for flat and steep-spectrum sources 
determined by ,Ricci et al., ( ,2004b) at 20 GHz, while the polari- 
sation angle is drawn randomly from a uniform distribution. 

Bright far-IR sources include those taken from the IRAS 
Point Source Catalog, with flux densities extrapolated in fre- 
quency using a grey-body spectrum, plus high-z galaxies 
strongly gravitationally lensed, based on the model by Perrotta 
et al. (2003). For fainter galaxies, which make up most of the 
cosmic infrared background, we adopted the model by Granato 
et al. (2004). Their clustering properties are modelled following 
Negrello et al. (2004). The implementation of their spatial dis- 
tribution is made using the algorithm by Gonzalez -Nuevo et al. 
(2005). The resulting power spectrum of CIB fluctuations turns 
out to be quite close to that measured by Planck Collaboration 



(Planck Collaboration XVIII 



201 1 1. The mean polarisation de- 



gree of individual infrared galaxies can be adjusted as an input 
parameter, the default value being 1%. The polarisation angle is 
randomly chosen from a uniform distribution. 

On the whole, these simulations should provide a useful tool 
for several purposes: to test data analysis pipelines for CMB ex- 
periments, identify the most convenient areas of the sky for spe- 
cific purposes (i.e., areas least contaminated by Galactic emis- 
sions), identify the optimal set of frequency channels for CMB 
temperature and polarisation experiments or for other studies 
(i.e., for studies of the cosmic infrared background), predict the 
levels of confusion noise, including the effect of clustering for a 
given frequency and angular resolution, and much more. Access 
to documented versions of the package and to reference simula- 
tions is available from a dedicated websit^ 
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